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I 


Introduction 


For  many  problems  associated  with  the  field  of  ballistic 
mechanics  experimental  methods  required  for  model  evaluation 
are  very  difficult  to  control  in  practice  and  are,  in  addition, 
very  expensive  to  perform.  Furthermore,  it  is  very  difficult, 
if  not  impossible,  to  obtain  time  dependent  information  experi¬ 
mentally  on  the  state  within  solid  materials.  Thus,  in  many 
instances  it  is  imperative  that  one  use  a  numerical  model  which 
simulates  the  physical  experiments  to  predict  the  required  in¬ 
formation. 

The  first  step  in  the  construction  of  the  model  is  to  define 
the  governing  differential  equations.  A  more  difficult  task  is 
to  choose  a  numeric nl  method  which  is  used  to  replace  the  differ¬ 
ential  equations  by  a  finite  difference  form;  the  resulting  equat¬ 
ions  must  be  solved  on  a  high  speed  digital  computer.  The  main 
features  that  one  strives  for  in  the  design  of  a  numerical  model 
are  the  attributes  of  accuracy,  economy  and  ease  of  operation  of 
the  code  for  the  designer-user.  High  velocity  impact  phenomena 
leads  to  very  severe  material  distortions.  The  computer  model 
presented  in  this  report  has  beer  used  to  predict  deformation 
states  for  problems  in  which  the  striking  velocity  lies  in  the 
range  of  0. 2-4.0  km/sec.  For  very  much  lower  impact  velocities 
computation  times  may  become  extreme;  there  does  not  appear  to  be 
a  restriction  on  the  method  for  higher  impact  velocities.  As  a 
result  of  difficulties  associated  with  Lagrange  methods  when  large 
distortions  are  present,  it  was  felt  that  an  Eulerian  formulation 
was  desirable. 

Eulerian  codes  are  characterized  by  a  mesh  which  is  fixed  in 
space  for  all  time;  the  material  "flows"  through  this  mesh.  As 
such  the  Eulerian  method  has  the  intrinsic  capability  of  represent¬ 
ing  numerical  solutions  to  problems  with  large  deformations  over 
long  periods  of  time  without  incurring  the  Lagrangian  penalty  of 
mesh  distortion.  As  a  result  of  maintaining  the  uniformity  of  the 
mesh  the  accuracy  of  the  method  is  preserved  and  the  need  of  opera¬ 
tor  intervention  for  rezoning  the  calculation  is  eliminated.  In 
this  paper  we  shall  describe  the  key  elements  of  an  Eulerian  code 
for  ballistic  problems  which  is  second  order  accurate. 

As  is  well  known,  and  unless  special  care  is  taken,  diffusion 
of  one  material  into  another  material  can  occur  in  an  Eulerian  for¬ 
mulation.  To  prevent  such  diffusion  between  materials  all  moving 
surfaces  which  bound  each  domain  are  defined  by  material  particles 
(called  tracer  particles)  that  are  taken  to  be  the  end  points  of 
piecewise  linear  segments  which  approximate  the  boundary.  These 
particles  can  move  freely  throughout  the  fixed  Eulerian  mesh;  their 
motion  is  determined  by  ordinary  differential  equations.  Although 
a  Lagrangian  calculation  is  necessarily  introduced,  the  usual  re¬ 
zoning  difficulties  associated  with  such  a  formulation  are  simply 


solved  since  the  problem  is  reduced  to  the  redistribution  of  mesh 
points  along  a  line  which  uses  arc  length  as  the  independent  vari¬ 
able.  On  this  string  one  need  not  conserve  quantities  when  re¬ 
zoning  is  performed;  indeed  rezoning  just  requires  maintaining 
uniform  spacing  of  the  particles.  At  free  surfaces  and  interfaces 
the  correct  boundary  conditions  are  imposed  with  no  integration 
performed  across  boundaries  separating  materials. 

The  asymptotic  accuracy  of  finite  difference  methods  is 
measured  by  their  order  of  accuracy.  For  a  first  order  scheme 
the  error  is  halved  when  the  mesh  spacing  is  halved.  For  a  sec¬ 
ond  order  method  halving  the  mesh  spacing  results  in  errors  that 
are  one  quarter  their  previous  value.  Numerous  comput  itional  ex¬ 
periments  have  verified  that  the  extra  accuracy  associated  with 
second  order  methods  more  than  compensates  for  the  additional  work 
required.  This  accuracy  can  be  utilized  in  two  ways.  One  can  fix 
the  mesh  spacings  based  on  computer  storage  facilities  and  specific 
problem  resolution  requirements;  the  higher  order  methods  will  yield 
greater  accuracy.  Alternatively,  one  can  fix  the  accuracy  desired. 
In  this  case  second  order  methods  will  need  fewer  zones  compared  to 
first  order  methods.  As  a  result  of  the  dependence  of  the  time  step 
on  the  mesh  spacing,  for  explicit  methods,  less  computer  time  is 
required  for  the  specified  accuracy. 

In  addition  to  being  second  order,  the  code  that  has  been  de¬ 
veloped,  called  SMITE,  uses  a  dissipative  scheme  in  divergence  form; 
hence  shock  waves  imbedded  within  the  flow  region  are  treated  auto¬ 
matically  and  with  the  correct  jump  conditions.  This  result  is  ob¬ 
tained  because  mass  momentum  and  energy  are  conserved  by  the  differ¬ 
ence  equations.  Our  model  has  been  formulated  in  the  cylindrical 
coordinates  (z,r). 
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II 


General  Modeling  Considerations 


The  differential  equations  describing  the  behavior  of  a 
continuous  media  involve  both  spatial  and  temporal  derivatives 
of  material  properties.  In  order  to  compute  the  rate  of  change 
of  certain  of  the  variables  one  requires  a  knowledge  of  certain 
spatial  derivatives. 

The  accuracy  of  the  temporal  changes  depends  not  only  on  the 
accuracy  of  the  difference  scheme  but  also  on  the  resolution  of 
the  mesh  in  the  neighborhood  of  the  regions  of  greatest  variation 
in  the  solution.  Such  regions  where  functions  can  vary  very  rapidly 
on  a  length  scale,  which  may  be  small  compared  to  a  characteristic 
length  in  the  region  where  the  differential  equations  are  solved, 
is  called  a  boundary  layer. 

It  is  obvious  that  imposing  a  fine  mesh  uniformly  over  the 
region  of  integration,  so  as  to  represent  the  function  numerically 
by  its  values  on  the  computational  mesh,  is  an  inefficient  proce¬ 
dure.  The  fine  grid  spacing  in  the  smooth  region  of  the  function 
is  unnecessarily  accurate  and  leads  to  large  computational  times. 

Let  r=a  be  a  line  in  a  region  where  the  function  is  slowly 
varying  while  r=0  be  a  line  in  the  boundary  layer.  Then  the  mesh 
spacing  at  r=0  should  be  less  than  the  spacing  at  r=a  by  a  factor 
d/a<l.  One  example  used  to  achieve  such  spacing  is  a  logarithmic 
variation  of  the  mesh  given  by 
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log  ( 


b+ (a-r ) 
b- (a-r) 
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,  ,b+a. 


(2.1) 


Here  L  is  a  measure  of  the  number  of  mesh  points  desired  in 
the  r  direction  and  d  is  the  length  scale  over  which  the  solution 
exhibits  its  largest  variation.  The  function  3(r)  is  the  change 
of  variables  such  that  6  is  monotonic  and  most  rapidly  increasing 
in  the  thin  boundary  layer  and  more  slowly  increasing,  with  in¬ 
creasing  r,  in  the  region  where  the  solution  is  smooth,  with  a 
suitable  choice  of  6,  such  as  the  example  (2.1),  there  will  no 
longer  appear  in  the  functions  which  are  approximated  on  the  differ¬ 
ence  mesh,  and  which  are  considered  to  depend  on  the  independent 
variable  0,  a  boundary  layer  type  of  structure. 

SMITE  has,  at  present,  a  transformation  for  z  which  introduces 
the  new  coordinate  a  through 

a  =  D  (z  -  g-h(r)]  ;  (2.2) 
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. . 


this  transformation  has  been  used  primarily  for  modeling  conical 
shaped  charges.  The  slope  cf  a  cone  in  the  physical  plane  is  q 
while  h  is  used  as  a  measure  of  the  size  of  the  cap  at  the  vertex 
of  the  cone.  In  all  subsequent  discussion  q  and  h  are  chosen  such 
that  a  =  Dz  is  the  transformation  with  D  a  conversion  for  centi¬ 
meters  into  coordinate  lines.  The  choice  of  parameters  used  in 
transformations  (2.1)  and  (2.2)  is  described  in  Section  VI. 

The  basic  equations  of  hydrodynamics  require,  at  each  time 
step,  the  evaluation  of  the  divergence  of  the  flux  of  the  con¬ 
served  physical  quantities,  i.e.  and  .  Under  the  condition 

of  coordinate  substitution  mentioned  above  in  which  3=6 (r),  I.e. 
Equation  (2.1),  the  divergence  of  the  flux  becomes 


3g  3q  33  ,  _  3f  3f  3a 

gr  =  37  =  #  37  '  fz  “  37  =  37  37 


We  find  that,  in  the  case  of  (2  1), 


33  _  dS 

37  “  77  = 


2  Lb 


which  shows  that  the  mesh  spacing  in  cylindrical  coordinates  is 
2  2 

proportional  to  b  -(a-r)  .  Hence,  as  r  varies  from  zero  to  a  the 

2  2  2 
mesh  spacing  varies  from  b  -a  ,  the  smallest  variation  to  b  which 

is  the  largest  variation 

The  introduction  of  such  coordinate  transformations  introduces 
the  concept  of  a  computational  plane  (coordinates  a-3)  in  which  the 
grid  is  uniform  but  whose  image  in  the  physical  plane  (coordinate 
z-r)  is  dense  in  the  boundary  layer  region  and  gets  progressively 
sparse  in  regions  where  gradients  in  the  solution  are  small.  This 
results  in  economical  use  of  the  computing  time. 

Through  application  of  the  chain  rule,  it  is  possible  to  trans¬ 
form  the  divergence  free  equations  from  the  physical  space  z-r  to 
the  computational  plane  a-3  in  such  a  manner  that  the  new  equations 
are  still  in  divergence  free  form.  A  direct  consequence  is  that 
internal  shocks  will  be  computed  with  the  correct  jump  conditions 
in  the  computational  plane. 

The  difference  schemes  that  we  use  to  approximate  the  diver¬ 
gence  form  of  the  differential  equations  requires  a  nine  point 
rectangular  lattice  on  the  initial  data  plane  at  time  level,  t. 
Symmetrically  placed  about  the  point  (i,j),  at  which  we  wish  to  find 
the  functions  which  constitute  the  solution,  are  eight  nearest 


neighbors  whose  position  is  defined  by  translations  about  (iAz,  jAr) 
by  ±Az  and  ±Ar.  Nine  point  stencils  which  are  entirely  inside  the 
domain  of  integration  at  time  t  are  updated  to  time  t+At  on  the 
basis  of  a  two  step  difference  algorithm  (see  Section  4). 

For  mesh  points  near  boundaries  the  nine  point  stencil  required 
by  the  two  step  method  is  no  longer  completely  contained  within  the 
domain  of  integration.  In  this  case  we  no  longer  consider  the  di¬ 
vergence  form  for  the  differential  equations  but  instead  solve  a 
quasilinear  system  of  equations  derived  from  the  conservation  equa¬ 
tions  for  the  unknown  functions.  The  solution  at  these  points  is 
updated  by  a  non-linear  version  of  the  Lax-Wendroff  method.  This 
method  allows  one  to  take  into  account  in  the  difference  equations 
the  non-centered  spacing  required  as  a  result  of  the  boundary  cross¬ 
ing  between  adjacent  points  of  the  stencil.  This  non-centering  occurs 
in  the  difference  formulas  since  data  along  the  boundary  must  now 
be  used  instead  of  data  at  regular  mesh  points.  As  any  of  the  eight 
nearest  neighbors  can  be  excised  from  the  stencil  by  the  boundary, 
there  exists  256  possible  ways  of  updating  the  solution  at  these 
"irregular  mesh  points".  Since  many  of  these  truncated  stencils  are 
mirror  images  of  each  other,  the  required  logic  can  be  minimized. 

In  order  to  calculate  the  non-centered  space  differences  to 
the  first  and  second  derivatives,  we  must  know  the  position  of  the 
boundary  together  with  the  value  of  the  dependent  variables  at  all 
crossings  of  the  boundary  with  coordinate  mesh  lines.  This  is 
accomplished  by  interpolation  from  the  material  particles  defining 
the  boundary,  care  being  taken  to  account  for  possible  multiple 
crossings  of  the  boundary  with  a  given  coordinate  line.  Once  all 
spatial  differences  are  known  to  second  order  the  solution  at  time 
t+At  is  calculated  by  a  truncated  Taylor  Series  for  the  unknown 
vector  function  W(t+At)  about  time  t,  keeping  all  terms  through 
second  order. 

The  moving  boundaries  consist  of  segments  which  satisfy  either 
free  surface  conditions  or  contact  discontinuity  conditions.  In 
addition,  there  may  exist  fixed  boundaries  some  of  which  are  lines 
of  symmetry;  here  we  use  the  usual  reflection  conditions  to  advance 
the  solution.  At  such  fixed  boundaries  the  solution  is  reflected 
across  the  boundary  to  image  points  which  allow  the  use  of  regular 
difference  equations  on  the  boundary.  All  surfaces,  which  move 
through  the  Eulerian  mesh,  are  marked  by  material  particles  which 
obtain  velocity  components  from  the  interior  of  the  domain;  their 
motion  is  governed  by  Lagrangian  equations  which  are  integrated  at 
each  time  step  to  predict  the  new  boundary  position. 

In  the  next  sections  we  make  a  more  detailed  discussion  of 
the  above  general  considerations. 
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Basic  Differential  Equations 


a)  Axisymmefric  Elastic  Model 

When  a  material  supports  shear  stresses,  it  is  necessary  to 
include,  in  addition  to  the  pressure  forces,  terms  which  account 
for  the  presence  of  these  stresses.  The  equations  of  motion  for 
such  a  material  can  be  derived  by  applying  the  physical  laws  de¬ 
scribing  the  conservation  of  mass,  momentum  and  energy  to  a  finite 
element  of  the  material  body.  In  addition,  a  statement  of  the 
stress-strain  relationship  of  the  material  is  required.  For  this 
paper  a  linear  theory  is  assumed,  i,e.  material  bodies  will  satisfy 
Hooke's  law.  Then  these  laws  may  be  usefully  written  as  a  set  of 
partial  differential  equations  in  a  cylindrical  coordinate  system, 
as  follows.  If  the  substantial  or  particle  derivative  is  defined 
as 


d  _  3 

at  "  5T 


u 


a 

Tz 


+  v 


3 

TE 


then  the  conservation  of  mass  can  be  written  in  terms  of  the  density 
(the  mass  per  unit  volume)  p,  and  the  divergence  of  the  velocity 
field,  with  components  u  and  v  in  the  z  and  r  direction  respectively, 
as 


dp 

dt 


3v 

3r 


+  1) 
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(3.a.l) 


The  two  momentum  laws 
components  x^.  The  axial 


reflect  the  appearance  of  the  stress 
momentum  equation  is 
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3r 


(3. a. 2) 


and  the  radial  momentum  equation  is 


dv  _  3t12  l  3t22  j  T22-T33 
pdt  3z  3r  r 


(3. a. 3) 


The  evolution  equation  for  the  internal  energy,  e,  per  unit 
volume  is  given  by 


de 

p3t 


lH  +  T  IX  +  T  (IX  +  3u. 

C11  3 z  T 22  ar  t12  l9z  ar' 


v(Tll+T22+3p) 
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(3. a. 4) 
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The  stresses  required  in  the  above  relations  must  be  obtained 
from  the  strains  and  strain  rates „  The  linear  stress-strain  laws, 
with  correction  for  rotation,  are  usually  written  in  terms  of  de¬ 
viator  stresses  S^j.  The  rate  of  change  of  the  stress  component 

Si;j  are  given  in  terms  of  the  strain  rate,  6^,  via 

-!£-?>  +  *i2  (3-a-5) 


dS12  ,3u  .  3v,  S11"S22  ,3u  3v. 
3t V  (57  +  37) - 2 (57  "  3z> 


(3. a. 6) 
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(3. a. 7) 


2u  ,3u  3v  2v. 

3  3  z  +  3r  r ' 


(3. a. 8) 


The  above  Hookian  laws.  Equations  (3 „a. 5) - (3.a. 8) ,  are  connected 
to  the  evolution  laws  (3.a.l)-{3„a.4)  by  the  algebraic  conditions 


T  .  .  *  S  .  .  -  p  6  .  . 
l]  l]  v  13 


6ij  “  1  for  1  =  j 

*  0  otherwise 


(3. a. 9-3. a. 11) 


The  pressure  p  is  related  to  the  density  p  and  specific  in¬ 
ternal  energy  e  through  the  equation  of  state 


p  =  P (p,e) 


(3. a. 12) 


The  above  set  form  a  system  of  twelve  equations  for  the  twelve  un¬ 
knowns  p,u,v,e,p,i11,T12,T22,T33,S11,S22  and  S33„ 

At  this  point  we  show  that  Equations  (3. a. 1)- (3. a. 12)  form  a 
system  which  is  not  self  consistent.  To  see  this  add  Equations 
(3. a. 5),  (3. a. 7)  and  (3. a. 8).  It  is  clear  then  that  the  sum 


Z  S..  satisfies 


dt  (S11  +  S22  +  S33)  “ 


0 


(3. a. 13) 


which  implies  that  the  sum  of  these  stress  deviators  is  a  constant 
of  the  motion  of  the  materia] ;  without  loss  of  generality  this 
constant  can  be  taken  to  be  ze»  t  for  at  t=0  each  S^=0.  Thus 


3 

I  S..  =  0  (3. a. 14) 

i=l  11 


for  all  time. 

Now  if  we  sum  Equations  (3. a. 9)  through  (3. a. 11)  for  i=j,  we 
obtain  the  relation 


T11  +  T 22  +  T  33  "  S11  +  S22  +  S33  ~  3p  (3,a,15) 
which  yields,  after  satisfying  Equation  (3. a. 14), 


T11  +  X22  +  C33  “  "  3p 


(3. a. 16) 


Equation  (3.a.l6)  states  that  the  pressure  p  is  determined  by  the 
mean  of  the  stress  tensor.  This  is  a  contradiction  of  Equation 
(3. a. 12)  which  states  the  pressure  is  a  function  only  of  the  density 
and  internal  energy. 

Hooke's  laws  can  be  written  in  the  form 


x .  .  =  2ue .  .  +  AZe .  . 
ii  ii  j  33 


i  =  1,2,3 


(3. a. 17) 


Tij  =  2ueij 


Here  p  is  the  shear  modulus  of  the  material, 
and  the  strain  e^  is  defined  by 


e 


ij 


1 

2 


‘SX- 


A  is  a  Lame  constant 


(3. a. 18) 


The  displacements  are  x^. 


Differentiating  Equations  (3. a. 17)  and  (3. a. 10)  with  respect 
to  time  yields 


Tii 


Tij 


2Meii  +  Xljejj 


2yeij 


(3. a. 19) 


The  corresponding  strain  rate  tensor  is  then  given  in  terms  of  the 
velocity  gradient, 


3u.  3u. 


®ij  “  7  (7x7  +  W 


with  and  Uj  the  components  of  the  velocity. 


If  Equation  (3. a. 16)  is  differentiated  with  respect  to  time 
we  can  compute  p;  using  Equation  (3. a. 19)  for  we  obtain: 

•  1  • 

p  -  -  i  E'ii 

■  -  j  [2p£iii +  3iESii] 

=  -  2^+-~  Ee.^  =  -  k  div  u 

-■  k  +  If  +  f>  (3-a-20) 


We  have  defined  the  constant  k,  the  bulk  modulus,  in  terms  of  the 
Lam£  constants  y  and  X. 

Clearly  Equation  (3. a. 20),  obtained  from  Hooke's  law  is  not 
consistent  with  Equation  (3. a. 12).  In  fact  Equation  (3. a. 20)  to¬ 
gether  with  Equation  (3.a.l)  yields 


dp  _  k  dp 

at "  p  3t 

which  can  be  integrated  from  the  lower  limit  (p,p)  =  (0,pQ)  to  the 
state  <p,p);  this  yields  for  the  pressure 

p  =  kLn (1+p)  ,  y  =  -  1  (3. a. 21) 

p0 


Wilkins1  points  out  that  Equation  (3. a. 21)  can  be  expanded  in  a 
Taylor  series  for  p/pg'tl  to  obtain  a  polynomial  for  p  as  a  function 

of  u.  Thus,  Equation  (3. a. 12)  should,  for  consistency,  be  close  to 
Equation  (3. a. 21)  but  obviously  cannot  in  general  be  the  same.  To 
be  consistent  then,  there  is  really  no  degree  of  freedom  in  the 
choice  of  the  form  for  an  equation  of  state  Equation  (3. a. 12)  in  the 
system  (3. a. 1)- (3. a. 12)  . 

One  possible  way  to  avoid  the  above  inconsistency  is  to  re¬ 
place  the  stresses  t ^ .  appearing  in  Equations  (3. a. 1)- (3. a. 4)  by 
the  deviatoric  stresses  and  the  pressure  p  through  the  use  of 

Equations  (3.a.9)-(3.a.ll)  and  then  eliminate  Equations  (3. a. 9)- 
(3. a. 11).  The  stress  deviator  S,..  can  also  be  eliminated  by  using 
Equation  (3. a. 14).  The  resultingJset  of  equations  can  be  written  as 


dp 

at  " ' 

3  u  3  v 

p (37  +  37 

+  -) 
r ' 

(3. a. 22) 

du 

3£  3  S11 

,  3S12  S12 

dt 

3  z  3  z 

dv 

dt  ” 

3£  3S12 

3r  3  z 

,  3S22  2S22+S11 

3r  r 

(3. a. 24) 

de 

at  = 

(Sll“p)  H 

+  (S22_P)  tr 

+ 

S12  <H + 

3  u,  v<Sll+S22+p) 

3r J  r 

(3. a. 25) 

dSn 

dt 

hi  (22“  - 

3  1  3Z 

n  _  v,  s  - 

3r  r'  &12  *3r 

3  z; 

(3. a. 26) 

dS12  . 
dt 

, 3u  ,  3v. 
p  3r  3Z5 

Sll-S22  ,3u  _  3v, 
2  l3r  3  z 

1 

(3. a. 27) 

dS22  _ 
dt 

hi  ( 2—  - 

3u  v,  „  ,3u  _ 

3 z  _  r'  b12  ' 3r 

— ) 

3  z ' 

(3. a. 28) 

P  =  P( 

p,e) 

(3. a. 29) 

We  now  have  eight  equations  for  the  eight  unknowns  p,u,v,e, 
S11'S12'S22  and  p>  Equations  (3. a. 22)- {3. a. 29)  are  consistent. 

The  only  disadvantage  of  this  set  of  eauations  is  that  Hooke's 
law  of  linear  elasticity  can  not  be  recovered.  This  is  due  to  the 
fact  that  the  pressure  in  Equation  (3. a. 29)  is  the  thermodynamic 
pressure  while  Hooke's  law  does  not  attempt  to  include  thermo¬ 
dynamic  effects.  If  nonadiabatic  terms  are  included  so  as  to  pro¬ 
duce  a  modified  Hooke's  law  then  the  relation  (3. a. 21)  would  be 
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modified  to  include  temperature  effects  while  the  stress  deviators 
Sij  would  be  unchanged.  Hence  the  quantities  in  Equations 

(3. a. 22)- (3. a. 29)  are  the  stress  deviators  and  an  additional  constant, 
i.e.  k  in  Equation  (3. a. 20)  (in  addition  to  the  shear  modulus  p) , 
usually  encountered  in  two  dimensional  linear  elasticity,  is  elim¬ 
inated  by  the  inclusion  of  an  equation  of  state  (3. a. 29). 

b)  A  Plastic  Model 


Equations  (3. a. 26)- (3. a. 28)  express  the  stress-strain  rela¬ 
tions  for  a  material  behaving  with  linear  elastic  properties.  Be¬ 
fore  we  proceed,  we  first  rewrite  these  equations  in  a  more  convenient 

form.  We  define  the  derivative  ^  to  be  a  tensor  operator  unaffected 

by  rotations.  Then,  we  construct  this  operator  from  the  substantial 
derivative  of  the  stress  deviators 


DS11 

dSll 

e  /3u  _ 

S12  (T?  T£> 

Dt 

3t 

DS12 

dS12  + 

Sll“S22  ,3u  3v, 

- 5 -  (37  “  TV 

Dt 

"  3t  + 

DS22 

dS22  4 

c  /3u  _  3v. 

S12  (Tr  TV 

Dt 

"  at—  + 

(3.b.l) 


The  time  derivatives  of  the  strain  deviators  can  be  written  using 
the  definitions: 


'll 


3u  1  ,3u  3v  .  v. 
7*  “  I  (3¥  +  37  +  r> 


2  3u  1  3v  1  v 


(3.b.2) 


•  _  3v  1  ,3u  3v  v. 
e22  3r  “  1  3z  +  3r  +  r' 


1  3u  2  3v  1  v 
T  3z  J  jr  "  T  r 


i  =  i  iSH.  + 
e12  2  lsr  3ZJ 


Then,  Equations  (3.a.26)-(3.a.28)  can  be  written  so  that  the 
deviator  stress  components  of  the  stress  tensor  are  obtained  from 
the  deviator  strain  components  of  the  strain  tensor,  i.e.. 
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(3.b.3) 


DS11 

Dt 

2Well 

DSi2 

Dt 

2vi12 

DS22 

• 

Dt - 

2pe22 

Equations  (3.b.3)  are  applicable  only  in  the  elastic  region  of  flow. 
In  general  a  material  which  exhibits  a  linear  variation  of  strain 
with  stress  is  called  linear  elastic.  However  at  the  proportional 
limit  the  strain  may  increase  more  rapidly  with  increasing  stress. 

In  this  region,  the  material  deforms  plastically.  If  the  strain  is 
allowed  to  increase  with  no  increase  in  the  stress  the  material  is 
called  perfectly  plastic.  If  some  variation  in  stress  occurs  the 
material  is  experiencing  work  hardening. 


A  material  when  exposed  to  external  loading  can  experience 
permanent  deformation  as  stresses  exceed  certain  characteristic 
limits  of  the  material.  A  tacit  assumption  is  made  in  elastic 
theory:  the  assumption  that  a  scalar  function  f (t^ j , e?j ,w) ,  called 

a  yield  function,  exists.  Arguments  T^j,  e?j  and  oj  correspond  to 

the  stress  state,  the  plastic  strain  and  a  measure  of  the  loading 
history  respectively. 


The  equation 


f  -  0 

represents  a  surface  in  stress  space;  for  f  <  0  the  change  in  plastic 
deformation  is  zero  while  only  when  f=0  is  plastic  deformation 
allowed  to  occur.  If  the  material  properties  are  independent  of 
strain  rate, f>0  has  no  meaning.  In  the  plastic  region,  in  place  of 
System  (3.b.3),  we  invoke  the  Prandtl-Reuss2  formulation  for  plastic 
flow. 

In  a  mixed  elastic  plastic  flow  material,  System  (3.b.3)  applies 
whenever 

Z  S?j  =  2  (S^  +  S^  +  sn  S22  +  S^2)  <  2K2  (3.b.4) 

1  Si,  j  <_3 

2 

with  K  a  constant  of  the  material.  However,  whenever  the  von  Mises 
yield  condition,  based  on  the  assumed  form  for  the  yield  function 
f=f(LS^j),  requires  that 

ZS2j  >  2K2  (3.b. 5) 

be  satisfied,  then  System  (3.b.3)  is  replaced  by  a  viscoelastic 
model  system  patented  after  a  viscoelastic  constitutive  relation 
of  the  form 
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(3,b.6) 


■  2|i 


Bt  “  (e*12  "  XS12) 


Bt  “  to  d 22  ~  XS22) 

Now  the  constant  X  is  determined  by  requiring  equality  in  the 
von  Mises  yield  criterion  (3.b.5)  rather  than  setting  it  to  1/2  pv , 
v  the  kinematic  viscosity  of  the  material.  Multiply  each  Equation 
(3.b.6)  by  S. .  and  sum: 


i  DS..  ■) 

7  ^BF-1  =  to  (ES^e*  -  XSS..) 


(3.b.7) 


Now  use  the  fact  that  Es2^  *  2K2  =  constant.  Equation  (3.b.7)  can 
then  be  solved  explicitly  for  X  , 

X  =  73  ZSije*ij  (3.b.8) 


In  cylindrical  coordinates  (3.b.8)  can  be  expressed  as 

i* 

X  =  — — <  s, ,  [u  -  1/3  (u  +  v  +  -) 

2K2  11  l_z  '  '  zz  rr  r'J 


[vr  -  1/3  <uzz  +  vrr  +  J>] 


+  S12  (ur  +  vz)  -  (S1]_  +  S22) 


|  -  -  1/3  (u  +  v  + 
l_r  '  '  zz  rr 


’  2K2  [Sl1  “ 


ur+vz 

zz  +  S22  vrr  +  S12  <— 1 1 


-  (sn  +  s22) 


(3.b.8* ) 


Here  we  have  used  the  notation  “  uz»  etc. 
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In  this  way  both  the  elastic  and  plastic  regions  can  be  de¬ 
scribed  by  Equation  (3.b.6).  The  prescription  is 


(rsL  <  2K2 


Elastic 


or  IS 


ij 


2K  and 


ESijEij<_0  (unloading) 


(3.b. 9a) 


Kj 


2K 


Plastic-/ 


A  = 


2K 


ESij*ij>0 


(3.b.9b) 


Equations  (3.b.9a)  and  (3.b.9b)  show  that,  in  the  plastic 
region,  if  one  begins  on  the  material  yield  surface,  and  in  the 
absence  of  unloading,  then  one  remains  on  the  yield  surface. 

Because  of  the  complicated  boundary  conditions  together  with  the 
nonlinearities  of  Equations  (3.b.9a)  and  (3.b.9b)  this  system  must 
be  solved  numerically.  If  one  uses  a  finite  difference  technique 
in  the  plastic  region  then  the  truncation  errors  inherent  in  any 
difference  scheme  will  result  in  a  set  of  deviator ic  stresses  which 
no  longer  lie  on  the  yield  surface.  It  is  therefore  necessary  to 
change  the  finite  difference  schemes  in  the  plastic  region  to  insure 
that  unloading  does  not  occur  due  to  truncation  errors.  Thus,  A 
in  the  numerical  method  will  not  strictly  be  determined  by  Equation 
(3.b.9b)  but  instead  the  derivation  of  this  formula  will  be  used  to 
force  the  yield  condition  to  be  satisfied  numerically. 

In  order  to  describe  the  method  used,  which  is  second  order 
accurate,  we  assume  that  a  solution  is  known  at  time  t  and  we  wish  to 
determine  the  solution  at  time  t+At.  The  solution  to  Equation  (3.b.l0) 
with  A=0  (i.e.  the  elastic  case)  will  be  denoted  by  s| .  .  Then,  by  using 
a  backward  Taylor  Series  in  time  one  has  -1 

Sij(t)  =  Sij<t+At>  “  At  Si;j(t+At) 

2 

+  -(-At)_  s.  .  (t+At)  +  01  ( (At)  3 ) 
or  J 

Sij(t+At)  =  Si;.(t)  +  At  (t+At) 

2 

-  sV  (t+At)  +  e/((At)3)  (3.b.l0) 
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Using  the  differential  Equation  (3.b.6)  in  Equation  (3.b.l0)  yields 
Sij(t+dt)  -  Si;1(t)  +  2Aty  [^(t+At) 

-  A(t+At)Si;i  (t+At)]  -  2y  ^  (t+At) 

A(t+At)SH,  <t+At)-A(t)S,  .  (t)-j 
- ^t-  - U - J  (3.b.ll) 

Or  introducing  the  elastic  deviatoric  stresses,  S? . ,  Equation  (3.b.ll) 
can  be  written  as  -1 

Si  j  (t+At)  «  (t+At)  -  2Atn  X  (t+At)  (t+At) 

+  Aty  (t+At) (t+At)-  X(t)Si;j  (t)|(3.b.l2) 

Because  all  terms  containing  S^j  (t+At)  are  linear,  we  may  solve 
directly  for  the  predicted  deviatoric  stress  at  the  advanced  time 
level  via 

Si j  (t+At)  =  ciJs^t+At)  -  AtpX(t)Si;j  (t)J  (3 „b.  13) 

All  terms  on  the  right  hand  side  of  Equation  (3.b.l3)  are  known 
except  for  a.  We  determine  a  by  requiring  the  S^j  (t+At)  to  lie  on 

the  yield  surface.  As  before  we  square  Equation  (3.b.l3)  and  sum 
over  i  and  j .  Then 

2K2  =  .1.  S2  .(t+At)  =  a2  E.  fsf  . (t+At) 

"i  2 

-  (At)yX  (t)sij  (tW 
Solving  for  a  we  obtain 


a 


K 


(t+At)-AtyX  (t) S± j (t) 


0 


2 


(3 .b. 14) 


To  sum  up,  the  procedure  for  solving  Equation  (3.b.6)  is 
given  by  the  following  three  step  algorithm: 
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i)  Determine  S®;(t+At)  by  solving  Equation  (3.b.6),  by  any 
second  order  method,  with  A«0. 

ii)  Test  if  £  Js® j  (t+At^  ^  <  2K^.  If  true,  set 

j (t+At)  =  si j (t+At)  otherwise  determine  a  from  Equation  (3.b.l4). 

iii)  Finally  solve  for  the  deviatoric  stresses  at  the  advanced 
time  level  using 

Si  j  (t+At)  =a|s®j(t+At)  -  (At)uX(t)SiJ  (t)j  (3.b.l5) 

with 

A  ’  <“’>  E 

1  lm  13 
1  <  n  <  3 


c)  Transformed  Differential  Equation 

The  partial  differential  equations  described  in  section  (3a), 
Equations  (3. a. 22)- (3. a. 28) ,  can  be  written  in  quasilinear  form,  i.e, 

w  +  A  w  +  B  w  +  i-  Cw  =  0  (3.C.1) 

u  z  r  r 


. -  -  ■ 


iAJtS'-Uti. 


i mi  ly^iha'ae.iniavtrivr 


A,  B  and  C  are  7x7  matriciesi  their  entries  are  given  below: 
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where  = 

Sij  ”  p^ij  and  ^nallV 

A 

0 

0 

0 

0 

0 

0 

/  0 

0 

0 

0 

0 

-1/P 

0 

0 

0 

0 

0 

-1/P 

0 

-2/p 

O  0 

Sll+S22+P 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

\  0 

0 

0 

0 

0 

0 

0 

\  0 

0 

0 

0 

0 

0 

In  Equation 

(3.C.1)  w  is 

the  seven 

vector 

I  p  \ 

I  u  \ 

V 


If  one  introduces  a  general  transformation 

a  =  a (z , r)  (3.C.3) 

g  =  e(z#r) 


we  can,  by  the  chain  rule  rewrite  (3.C.1)  in  the  oi-g  plane: 

w.  +  (Aa  +  Ba  )  w  +  (Bg  +  Ag  )  w 
z  z  r  a  r  z  p 

*  ?(«.6)CW-°-  (3’c-41 
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r  ,ri^  rt 


In  order  to  solve  the  system  (3.c.4)  by  the  second  order  Lax- 

Wendrof £"  method,  which  uses  s  Taylor  expansion  for  the  solution 
vector  w(t+At)  about  the  initial  data  w(t)  via 

2 

w(t+At)  =  w (t)  +  it  wt  +  wfct  +  ©'(A t3)  (3.C.5) 

It  is  necessary  to  compute  the  second  time  derivative  w.  ;  the 
first  time  derivative  is  obtained  directly  from  (3.C.4): 

It  is  sufficient  to  show  how  this  is  accomplished  for  a  single 
component  of  w.  For  example  the  equation  for  the  density,  from 
Equation  (3.C.1),  is  just 

P,.  +  UP„  +  vP_  +  Pu„  +  Pv,  +  ^  =  0  (3.C.6) 

t  z  r  z  r  r 

The  counterpart  of  Equation  (3.C.6),  in  the  (a-B)  plane  is 

just 

Pt  +  P(a2ua  +  arva  +  Brve  +  ezuB)  +  (u«z  +  Var)Pa 
+  (v6r  +  uBz)pb  +  *  0  (3.C.7) 


Under  the  assumption  that  the  coordinate  system  is  independent  of 
time  we  may  compute  the  second  time  derivative  of  the  density  from 
Equation  (3.C.7): 


'tt 


p.  (a  u  +  a  v 
t  z  a  r  a 


SrvB  + 


-  p(azuta  +  Vta  + 


BrvtB  + 


SzutB) 


Pta(uaz  +  var)  (3.C.8) 

(v3r  +  uBz)Pt6 
=  0 

Now  the  terms  on  the  right  hand  side  of  Equation  (3.C.8)  can  be 
obtained  from  Equation  (3.C.4)  by  taking  the  appropriate  space 
differential.  For  example  the  terms  involving  pfca  (and  pfcB>  are 

evaluated  by  operating  on  the  first  (second)  component  of  System 


“  (vt8r  +  ut6z)p8  “ 
-  ?(a,B)<ptv  +  VPt> 
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(3.C.4)  with  the  differential  operators  ^  (y^-)  •  Cross  deriva¬ 
tives  of  the  velocity  components  are  evaluated  by  taking  appropri¬ 
ate  space  derivatives  of  the  second  and  third  components  of  the 
System  (3.C.4).  In  this  way  we  find  that  only  spatial  derivatives 
appear  in  the  differential  relation  describing  the  second  time  deri¬ 
vative  of  the  density. 


3  3  3  ^  3^  g  2 

ptt  =  P(^  '  n  '  ^2  '  TcuTB-  ' 


(3.C.9) 


Here  P  is  a  second  order  differential  operator  acting  on  w. 

In  this  manner  all  components  of  the  vector  wfct  may  be  computed. 

In  the  present  study  we  considered  the  class  of  transformations 
to  be  restricted  under  the  assumptions 


az  =  1 


3z  =  0 


(3.C.10) 


Thus  3  is  associated  only  with  r  while  a  can  depend  on  both  r  and  z, 
Equation  (3.C.7)  is  then  simplified: 


p.  =  -  (u  +  a  v)p  -  P  (u  +  a  v  ) 
Kt  r  a  'a  r  or 


(3.C.11) 


-  Sr(Pve  +  vp6) 


pv 

r (a, 3) 


We  now  write  down  the  remaining  components  of  System  (3.C.4) 
under  condition  (3.C.10): 


The  axial  momentum  equation  is  used  to  obtain  the  rate  of 
change  of  the  axial  velocity  through 

u.  =  -  (u  +  a_v)  u  +  ^-(x,.  +arSn  ) 
t  r  a  p  11 ,  a  r  12 ,  a 


^12  8  ^12 
-  sr  (“'-a  -  +  — 


(3.C.12) 


The  radial  momentum  equation  is  used  to  oL  iin  the  rate  of 
change  of  the  radial  velocity  via 


v .  =  -  (u  +  a  v)v  +  -{S,.,  +  a  x,,  ) 

t  r  a  p!2,a  x  22,  a 


3r(v-ve 


-v„  -  + 


2S22+S11 


(3.C.13) 
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The  rate  of  change  of  the  specific  internal  energy  is  deter¬ 
mined  from 


'a  +  p  [(T11 

+  Vl2)ua 

+  V22)va] 

-  0r[ve0  (3.C.14) 

n 

v(p+Su+S22) 

!v6  +  S12ue2J 

pr 

The  final  three  equations  (3.Col5),  (3.C.16)  and  (3.C.17) 
determine  the  rate  of  change  of  the  deviatoric  stresses  S^,  S^2 

and  S22  respectively.  The  notation  corresponds  to  Szz,  S12 

to  Szr  and  S22  to  S^t 


Sll,t  (u  +  aj:v)Sll,a  +  (3”  +  arS12)ua 

(S12  J  ar  va  '  Br  [vSll,0  (3.C.15) 

.  c  u  +  2*i  v  1  -  2P  v 
S12US  3  @J  3  r 


>12,  t  -  -  <“  *  V"Sl2,a  *  »r 


sll"S22 

+  (  -V  -  + 


S11_S22  \ 
/ 


y)vu  -  6r  (vs12f| 


(3.C.16) 


’22,  t 


-  (u  +  arv)S22,c  ■ 
+  <S12  *  «r  3H'v. 


(|i  +  zrSu)ua 


(3.C.17) 


er(vS22g  +  S12u6  “  3^)  “  3^  r 


Equation  (3.C.12)  through  (3.C.17)  are  valid  for  6/  0.  For 
the  special  case  of  flow  on  the  axis  of  symmetry,  6=  0,  L'hospitals 
rule  when  applied  to  the  above  system.  Equations  (3.C.11)- (3.C.17)  , 
yields 


*  — 

t 

°uo  •  upa  •  2Br  PVB 

(3.c.lla) 

‘t 

+  5*11, a  +  “c 

(3.c.l2a) 

t  - 0 

t 

uea  +  ?  (Tllua  +  2f3rT22  V 

(3.c.l3a) 

(3.c.  14a) 

11,  t  = 

-  uSll,a  +  r  (ua  -  3rV 

(3.c.l5a) 

12, t  * 

0 

(3.c.l6a) 

22, t  = 

‘  uS22,a  "  T-K  -  PrV 

(3.c.l7a) 

System  (3.C.11)  through  (3.C.17)  together  with  (3.c.lla)  through 
(3.c.l7a)  are_the  differential  equations  solved  in  a  strip  of 
thickness  6^/2a,  A  the  spatial  step  size  near  the  boundary  of  the 
domain.  This  is  shown  in  the  figure  below  by  the  crosshatched  area. 


interpolation  annulus 


The  thin  annulus  directly  adjacent  to  the  boundary  is  a  region  where 
mesh  points  lie  too  close  to  the  boundary.  Since  the  stability  of 
the  finite  difference  solution  would  not  be  satisfied  in  this  region 
interpolation  between  the  boundary  data  and  interior  data  is  used 
to  update  the  solution. 

It  now  is  appropriate  to  describe  the  form  of  the  differential 
equations  used  for  the  solution  interior  to  the  domain  (the  region  I 
in  the  above  figure) .  Here  we  write  the  first  four  differential 
equations  in  conservation  form  choosing  the  entries  of  the  vector  w 
to  be  the  quantities  conserved  across  discontinuous  transitions,  i.e 

P 
pu 
pv 
£ 


w 


Here  E  is  the  sum  of  the  specific  internal  and  kinetic  energy,  the 
total  energy,  E  ■  P  (e  +  1/2 (u2  +  v2)). 

The  continuity  equation  is 

Pt  +  <pu)2  +  (pv)r  +  £1  =  o  (3.C.18) 

The  axial  momentum  equation  is 

(pu)t  +  (pu2  -  Tn)z  +  (puv  -  S12)r 

(puv-S. ,) 

+ - g  ±z  =  0  (3.C.19) 


The  radial  momentum  equation  is 


(pv)t  +  (puv  -  S12)z  +  (pv^  -  T22)r 


(pv  i ) 

+  - _£i — =  o  (3.C.20) 


Conservation  of  energy  requires  that  E  satisfy 

Bt  *  [<*  '  TU)»  -  Suv]2  *  -  T22)V-S12u] 


(E-t,,)v-S,,u 
+  - — - —  =  0 


r 

(3.C.21) 


The  stress  strain  relationships  are  not  relationships  that 
express  a  conservation  principle.  Hence  they  are  rewritten  here 
in  their  quasilinear  component  form: 


Sll,t  +  uSll,z 


3l2,t  +  uS12, z 


+  vSll,r  +  f^-2uz  +  vr  +  F 


+  S12  (Vur> 


(3.C.22) 


+  vS 


12 


,r  “  ^(ui  +  vz> 


SH"S22 


(u  -v J 


’22, t 


+  uS 


22, z 


+  vS 


1  -r  z 
2 


(3.C.23) 


22 


+  (-2v  +  u  +  -) 

,r  3  r  z  r' 


+  S12  (Vvz>  = 


0 


(3.C.24) 


The  specification  of  the  equation  of  state,  Equation  (3. a. 29), 
completes  the  system. 


Along  the  axis  of  symmetry  System  (3.C.18)  through  (3.C.24) 
are  redefined  by  the  application  of  L' hospital's  rule.  Hence  in 
z-r  coordinate  along  the  line  r«0,  we  obtain  the  system 


pt  +  <pu>z  +  2<pv>r  "  0  (3.c.l8a) 

(pu)  t  +  (pu2  -  T11)z  +  2  (Puv  -  S12)r  =  0  (3.c.l9a) 

(pv) t  =  0  ,  S^  +  2S22  =  0  (3.c.20a) 

Et  *  [(E  -  tu)U  -  S12v]j  ♦  2  [(E  -  t„)v 

-  S12uJ  r  =  0  (3.c. 21a) 

Sll,t  +  uSll,z  +  IP(vr  *  V  =  0  (3.c. 22a) 


S12,t  =  0 


(3.c.23a) 


S22 , t  +  uS22,z  +  Tp(uz  "  vr)  “  0 


(3.c. 24a) 


In  a  situation  where  the  transformation  to  the  a-g  plane  is 
the  identity  transformation,  the  System  (3.C.18)  through  (3.C.24) 
augmented  on  the  line  6=  r  =  0  with  System  (3.c.l8a)  through  (3.c.24a) 
would  be  the  complete  set  of  interior  equations  to  be  solved. 

However  for  transformations  used  in  the  present  work  (Equations 
(3.C.3)  subject  to  (3.C.10))  the  conservation  form  for  the  equations 
in  the  a- 8  plane  become 

(rBp)t  +  [Vpu  +  arpv>]  a  +  (pv)g  +  FIT  =  0  (3.C.25) 


(rgpu) 


«♦  0 


r^  (pu  -  +  ar(puv  -  S12) 


+  (puv  -  S, 


puv-S12 


)PUV-D 

8  +  ~ r6^ 


(3.C.26) 


(r8pV)t  +  [r6(pUV  "  S12  +  ar{pv2  '  T22)3  a 

I  2  \  ,  pv2”2S22”SH  _  „ 

+  r  -  T22j  6  *  —rr - ° 


pv 

+  - - —  =  0  (3.C.27) 
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The  above  four  equations  are  augmented  by  the  evolution  equations 
for  the  deviator ic  stresses  in  the  a- 8  plane,  Equations  (3.C.15), 
(3.C.16)  and  (3.C.17). 

Along  the  axis  of  symmetry,  6=0,  the  above  conservation  laws 
reduce  to 


(r6P)t  +  (r6pu)a  +  2  (pv|  8=  0 

(3.c. 25a) 

(repu)t  +  [r8(pu2  -  Tllj]a 

+  2  (puv  -  S12)6=  0 

(3.c. 26a) 

(repv)t  =  0  ,  2S22  +  8U  -  0 

(3.c. 27a) 

(r8E)t  +  [rB«E  '  T11)U|]  a  +  2  [(E  -  T22>  v 

-  S12uJ  e  =  0 

(3.c. 28a) 

The  deviatoric  stress  components  comp] eting  System  (3. c. 25a)- (3. c. 28a) 
are  determined  by  (3.c.l5a),  (3.c.l6a)  and  (3.c.l7a). 
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IV. 


Finite  Difference  Equations 


The  partial  differential  equations  described  in  Section  III 
have  been  written  in  both  conservation  form  and  quasilinear  form 
in  the  computational  (a-g)  plane.  This  is  necessary  since  one  set, 
the  conservation  form,  is  used  at  all  points  which  are  interior 
to  the  domain  and  contain  all  their  eight  nearest  neighbors  inter¬ 
ior  to  the  domain.  The  quasilinear  form  is  used  to  construct  the 
difference  scheme  to  be  used  at  interior  mesh  points  which  have  at 
least  one  of  the  nearest  neighbors  exterior  to  the  domain.  We  start 
the  discussion  with  the  main  difference  scheme  used,  the  two  step 
method  for  the  conservation  form  of  the  defining  partial  differen¬ 
tial  equations. 


a)  Two  Step  Method 

We  wish  to  solve  the  set  of  equations  defined  in  the  previous 
section,  Equations  (3.c. 18) - (3 .c . 24 ) ,  on  a  set  of  mesh  points 


a -  =  iAa  ,  i  =  0,1, . . . ,  I 


=  j A3  ,  j  =  0,1, ... ,  J 


(4.a.l) 


tn  =  nAt  ,  n  =  0,1, .. . 


For  convenience  we  introduce  the  notation  for  the  first  four 
equations  (3. c. 18-3. c. 21) .  Let  f,g  and  h  be  four  vectors  defined  by 


f  = 


and 


\ 


Pu 


pu  -i 


11 


puv-S 


12 


(E“Tll)u"S12V/ 


g  = 


/ 


ov 
PUV-S 


\ 


pv2-t 


12 


22 


(E-T22)v-Si2u 


h  = 


pv  "2S22_S11 
(E-t22) v_s12u 


/ 


(4„a. 2) 
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Then  with  wT  ■  (p,pu,pv,  E)  we  have 

wfc  +  fz  +  gr  +  i  h  -  0  r  +  0  (4. a. 3) 

wt  +  fz  +  (g+h)r  =0  r  =  0  (4. a. 3a) 


In  the  a-6  plane  (4. a. 3)  becomes/  after  application  of  the  chain 
rule 


*  «r  *.  *  6  A  +  ^ 


(F’  +  (F^)  +  g6  +  rS“ 

r  t  r  a  r 


(4. a. 4) 


'V 

where  f  =  f  +  a^g. 

Comparison  of  Equation  (4. a. 4)  with  the  component  forms. 
Equations  (3.C.25)- (3.C.28) ,  give  the  individual  entries  for  the 

'X/ 

new  flux  vector  f  and  g.  For  the  remainder  of  the  discussion  we 
drop  the  tilda  on  f. 

The  approximate  solution  is  called  V; 


V(ai'6j'tn)  =  Vij  %  w(ai'ej'tn)  (4.a.5) 

The  approximate  solution  is  written  as  a  two  step  difference 
equation.  Predicted  values  V  are  first  obtained  at  the  midpoints 
(ai+1y2,0 j+1y2't+nAt)  of  the  mesh  by  a  first  order  difference 

approximation.  These  values  are  then  used  to  obtain  a  second  order 
accurate  solution  at  regular  mesh  points.  Letting  \^  =  At/Aa 

and  A 2  =  At/A 6  the  finite  difference  equation  for  the  first  step 
is 


viJi/2,  i+i/a  ■  <v?zi,i+i *  *j+i,i  +  »;;i+i +  <)> 

-  i'»i<fKi.i+i  -  £2;i+i  *  <1:1,1-  <i> 


-  v>i2«slii,  i+i  -  ®Kl.i +  <i+i-  <i» 

-  *  <}+i 


+  <j'^'Bj+l/2> 


(4. a. 6a) 


Introducing  the  notation  ?-f(V)  the  second  step  is  defined  by 
the  finite  difference  equations 


•>n  ,,n-l  i  /  j  \  /  I  1  ,  "^n 

i,j  “  i,j  ‘  1/4  l(fi+l,j  '  fi-l,j  ?i+l/2, j+1/2 


-n-l 


n-1 


-  7n  +  Tn  -  Tn  ) 

1-1/2 , j+1/2  +  ri+l/2, j-1/2  ri”l/2, j-1/2' 


-i/*\  ,  n-1  „n-l  .  -n 

-  l/4X2(gi>j+1  -  +  ^i+i/2, j+1/2 


-n 


-n 


-n 


gi+l/2, j-1/2  +  gi-l/2, j+1/2  “  gi-l/2, j-1/2 


(4. a. 6b) 


-  i/44t(hj;J+i  +  hj;£.x  +  i/2(h;tl/Jfj+1/J 


+  h" 
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Stability  of  the  above  difference  scheme  is  achieved  if  an 
artificial  viscosity  Q  is  added  to  the  right  hand  side  of 
Equation  (4. a. 6b) 

0  -  K  |k1[luln>j  -  (V1+lij  -  vlf)>  -  lui  -  (VlfJ 

'  Vl-1. i D  *  X2  DVj+l  ‘  'i.  j'  -  Vi.  j> 

"  |vi,j  "  vi,j-ll  <Vi,j  "  Vi,j-1)]  f  (4. a. 7) 

'll  'U  (  •  — ■ 

where  u  and  v  are  equal  to  a  and  6  respectively  (Equation  (5.a.l(). 

The  time  step  At  is  kept  at  approximately  2/3  of  the  maximum  allow¬ 
able  CFL  value,  i.e..  At  *  .65  Atri^  .  We  compute  the  Courant- 
Friedrichs-Lewey  time  step  by  finding  n» 

n  =  max  j(u  +  +  1  C)/Az,  (v  +  BrC)/Ar| 

over  all  mesh  points.  The  maximum  time  step  is  then  AtCFL=  1/Dn. 

Equations  (3.C.22),  (3.C.23)  and  (3.C.24)  for  the  deviatoric 
stress  components  are  solved  in  an  entirely  analogous  manner. 

Now  let  w  —  (S^j  S ^ 2 *  ®22^y  ®  *  (P»u»v,e,S^2,  ^12f  ^22^  ' 


f  = 


I  4  \ 

JPU  ' 


-  liV 


\H 


( iuv\ 


-  yu 


/  jMV 


h  = 


(4. a. 8) 


^  /  \ 


and  define  the  matricies  a  and  b  as 


/ 


a  = 


0  - 


12 


Sll"S22 


(4. a. 9) 


V 


-S 


12 


■/ 


0 


0 


0 


0 


and 


-S 


12 


Sll"S22 


’12 


\ 


v  0  0 


0  v  0 


0  0  v 


Then  introduce  two  transform  matrices.  A  and  B  given  in  terms  of 
a  and  b,  A  =  a  +  (a+arb)/Br  and  B=b.  In  vector  notation  the  form 

used  to  generate  the  difference  scheme,  in  terms  of  the  above 
matrices,  with  £  =  (f  +  a^g)/^^  is  just 

(^)t*  £a  +  »6  *  A”a  +  %  *  ?6r  ’  0  <)-»-10> 

If  the  terms  with  coefficients  in  Equation  (4„a.l0)  when  put 
into  difference  form  are  centered,  the  same  two  step  algorithm 
(4. a. 6a)  and  (4. a. 6b)  results  for  the  stress  deviators  w. 


b)  One  Step  Method 

The  basis  for  the  one  step  algorithm  is  the  method  proposed  by 
Lax  and  Wendroff.  As  stated  in  Section  3,  a  Taylor  series  is  used 
to  determine  the  solution  at  time  t+At  from  known  initial  data  at 
time  t  via 


w  (t+At)  =  w  ( t)  +  At  wfc  +  wfct  +  er  (At3)  (4  .b.  1) 

In  the  last  section  we  described  how  the  system  of  equations  (3.C.11) 
(3.C.17)  is  solved  using  (4,b„l).  Reiterating  briefly  -  one  may 
solve  the  system  (4.b.l)  because  w,  is  determined  directly  from  the 
differential  equations.  In  order  to  evaluate  the  second  time  deri¬ 
vative  one  differentiates  these  same  equations  with  respect  to  time. 
Two  new  cross  derivatives  appear  -  namely  w  and  wg  ;  these  latter 
two  are  determined  by  directly  differentiating  the  same  system 
(3. c. 11)- (3. c. 17)  first  with  respect  to  a  and  then  differentiating 
with  respect  to  g.  Back  substitution  of  these  two  cross  derivatives 
into  the  wtt  equations  leads  to  a  differential  form  for  (4.b.l) 

the  right  hand  side  of  which  is  completely  independent  of  deriva¬ 
tives  with  respect  to  time. 

We  now  plan  to  be  quite  specific.  We  will  carry  out  all  the 
indicated  differentiation  in  the  a-B  plane  for  Equation  (4.b.l). 
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The  time  derivative  of  the  continuity  equation  (3.C.11)  is 

ptt  "  "  (u  +  arv)pat  “  <ut  +  arvt)pa  '  pt(ua  +  "rV 

"  p(uat  +  Vat*  ’  6r  (ptvB  +  pv6t  +  vpBt  (4*b‘2) 
+  vtpej-  i  (ptv  +  pvfc) 


Observe  the  appearance  of  at  and  6t  cross  derivatives.  They  will 
be  defined  shortly. 


The  time  derivative  of  the  momentum  equations  (3.C.12)  and 
(3.C.13)  are 

utt  =  “  (u  +  arv)uat  “  (ut  +  Vt)ua 


and 


p  [Tll,at  +  “rS12,at  "  pt(Tll,a  +  arS12,aj] 
0r  |vtufl  +  vu0„  -  -12*-3-  +  ^  S 


B  pt 


+  pr  {S12,t  “  p”  S12) 


p  ^2  12,6 


vtt  =  '  (u  +  arv)vat  -  (ut  +  “rvt)va  +  p  [S12, 

+  arT22,at  “  p^(S12,a  +  arx22,a)] 

-  6r  (  vt 


L_V  n  +  W  f 


.^t  +  £tT> 


'6  '  ”6t  p  '  p2  22 

*  f?<2S22,t  *  sU,t  ‘  FI2S22  *  Sll,> 


,b) 


(4 .b. 3) 


at 


(4 .b. 4) 


The  cross  derivatives  in  space  and  time  also  appear  in  the  above 
two  relations. 


The  second  time  derivative  of  the  internal  energy  is  the 
relative  complicated  expression 


ett  "  "  <u  +  arv)eat  "  <ut  +  arvt)ea 

+  ?  [<Tll,t  +  arS12,t)ua  +  <T11  +  arS 


12>uat 


+  (S12,t  +  arT22,t)vcx  +  (S12  +  “rT22)vat 


Kt 

p 

«(,u 

+ 

arS12)ua 

+  (S^2  +  a 

rT22)v< 

3r 

[vtee 

+ 

veet  "  £ 

(T22,tV6  + 

T22v6t 

si2,tue 

+ 

si2uet  - 

Pt 

P  (T22vB  + 

S12UB 

i 

P  r 

(p 

+  SU  +  ; 

S22  +  V(pt 

+  Sll, 

ft 

P 

v(p 

+  £ 

311  +  S22 

>] 

(4 . b. 5) 


The  term  pfc  is  evaluated  by  considering  p=p(t)  and  e=e(t)  in 
Equation  (3. a. 29);  hence  pfc  =  ppPt  +  peet* 

The  second  time  derivative  of  the  deviatoric  stress  is 

Sll,tt  =  -  (u  +  arv)Sll,at  "  (ut  +  arvt)Sll,a 

+  (  +  arS12)uat  +  arS12ftUa 

-  (S12  +  ar  -  S12>tva  I 

“  6r(VtSll,e  +  vSll,Bt  '  S12,tUB 
’  S12uet  +  ¥  V6t)  "  ¥  Vt 


(4 .  b.  6) 
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S12,tt  -  "  <u  +  arv)S12,at  "  <ut  +  Vt)S12,a 

+  (2P  -  Su  +  S22)uat  -  -§  (Sn<t 

'  S22,t)ua  +  !/2  (2P  +  Su  -  S22)vat 

+  1/2  (Sll,t  "  S22,t)va  'er  [vtS12,0 

+  vS12,0t  "  1/2  <2p  ~  sn  +  s22)u0t 
+  1/2  (Sn^t  -  S22,t)ueJ 


C4.  b.  7) 


S22,tt  -  *  <“  *  arv*S22,at  '  <“t  +  Vt»s22,a 

'  l¥  +  arS12)uat  -  “rS12,tuu  +  <S12 
+0lr  r)Vat  +  S12,tva  "  [vtS22,0+  vS22, 

+  S12,tu0  +  S12u0t  "  T-  v0t]  ~  T"4 


(4.b.8) 


On  the  axis  of  symmetry,  0=0,  Equations  (4.b.2)- (4.b.8) 
become 

ptt  =  -  upat  “  utPa  "  pt  ua  '  PV 


‘  2er<ptve  +  pv6t) 


U*-f  =  -  +  T  (T 


Pt 

-  — tn 


at  “to  p  '  ll,at  p  11, a 

2er  pt 

+  p  (S12,0t  "  ~  S12,0) 


(4.b.2a) 


(4.b.3a) 


Vtt  =  0 


(4.b.4a) 
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ett  “  -  ueat  -  utea  +  f  <Tll,tua  +  T11  uat 
pt  20r 

_  ~  TllUa)  +  T"  (T22,tV0  +  T22VSt  (4.b.5a) 

Pt 

p”  T22v6) 

511, tt  =  "  uSll,at  ~  utSll, a  +  T"  (uat  "  3rv6t)  (4*b*6a) 

512, tt  =  ®  (4.b.7a) 

S22,tt  =  "  uS22,at  ”  ut  S22,a  “  3^  (uat  “  8rv0t*  (4,b,8a) 


We  must  now  compute  the  at  and  St  cross  derivatives  which 
appear  in  the  right  hand  side  of  the  Systems  (4.b.2)- (4.b.8)  and 
(4.b. 2a) - (4 .b. 8a) ;  first  the  at  cross  derivatives  are  computed 
followed  by  the  St  derivatives. 

The  continuity  equation  yields  for  the  cross  derivative  of  the 
density  p^,  for  g>0. 


p  .  =  -  (u+o.  V)  p  -  2  (u  +  a_v  )p  -  p  (u  +  a  v  ) 
*at  r  Maa  a  r  a'  a  H  aa  r  aa 


VS  +  pVaS  +  Vap0  +  V 


Pas] 


(4.b.9) 


pv  +  p  v 
a  a 


while  on  the  axis  of  symmetry,  0=0,  the  density  satisfies 

p  .  =  -  up  -  2u  p  -  pu  -  2 0  ( p  v.  +  pv  )  (4.b.9a) 
'at  paa  a ’a  aa  r  'a  0  a0 


The  cross  derivative  of  the  axial  velocity  u  t  for  0>O,  satisfies, 


■  a  v)u 
r  aa 

(ua  + 

arV 

u 

a 

. .  +  a_S 

11, aa  r 

12,  aa 

-  fa 

P 

(Tll,a  + 

VaU0  +  VUa0 

''Ih*  -  fl2p 

-  fi2 

+ 

S12,0pa 

J 

P 

2 

P 
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HtffltfBtA  <Ht  iBiMMi  i 


while  for  8-0  (4.b.l0)  becomes 


-  uu  -  u  u  +  i  (x,. 

aa  a  a  p  11, aa 

20  P„ 

+  ~  (S12.a8  "  p“  S12,6) 


-  TT  T„  J 


(4 .b. 10a) 


The  cross  derivative  of  the  radial  velocity  va  for  8>o,  satisfies, 
vat  “  “  <u  +  arv)vaa  “  +  arva>  va 


—  Ts  +  a  t  -  . a  (c 

P  L*12,aa  T  r  22,00  p  *B12,a 

arT  22, a  j]  '  3r  (vav8  + 

T22,BPa\  ,  ,2S22,a+Sll,a 

“72  I  1  pr 


<2S22+Sll)Pa, 


(4 .b. 11) 


while  for  8=0  (4.b„ll)  becomes 


v  .  =  0 
at 


(4 .b. 11a) 


The  cross  derivative  of  the  internal  energy  eat  for  8>0,  satisfies, 


(u  +  a  v)e  -  (u  +  a  v  )e^ 
r  '  aa  a  r  a  a 


+  p  [(Tll,a  +  arS12,a)ua  +  (T11  +  “rS12)uo 
+  (S12,a  +  arT22,a)va  +  (S12  +  V22)vaa 


-  f  ((TH  +  arS12)ua  +  <S12  + 


arT22)va)]  (4*b*12) 


!r  [v«e8  +  Vea8  “  £(T22,av0  +  T22va0  +  S12,auB 


pa  “1  lT  (P+S11+S22) 

S12ua8J  “  “  (x22vB  +  S12U6))J  “  rLVa  P 

(pa  +Sll,a+S22,a)  Vpa  (P+S11+S22>1 

P  2  J 


(Pg  +Sll,g+S22,a)  vpa (P+Sll+S22>' 
P  P2 


+  v 


F 


while  for  0-0  ( 4 . b . 1 2 )  becomes 


-  ue  -ue  +  —  (t  , ,  u  +  t,,u 

ua  a  a  p  '  ll,a  a  11  aa 

P  26 

-  p2-TllUa)  +  ~r  <T22,avB 


(4 .b. 12a) 


+  T22va6  "  “  T22v6J 


Again,  the  pressure  derivative  pa  encountered  above  is  treated  in 
precisely  the  same  tashion  as  the  derivative  p.  described  earlier. 
The  cross  derivative  of  the  deviatoric  stress  S,,  for  B>0, 


satisfies, 


'11, at 


Sll,at  =  -  (u  +  arv)Sll,aa  “  (ua  +  arva)Sll,a 

4u 

+  (y-  +  «rS12)uaa  +  ar®12,aua 
‘  <S12  +  “r  I^aS  S12,avc*  “  Br(vaSll,f 


+  vSll,aB  “  S12,aUB  "  S12UaB  +  T~  va&) 


(4 ,b. 13) 


n  v 

in  _« 

3  r 


while  for  B=0  (4.b.l3)  becomes 


’ll, at 


uSll ,aa  “  uaSll,a  +  3  (uaa  "  0rvaB)  (4.b.l3a) 


The  cross  derivative  of  the  deviatoric  stress  S. -  for  B>0, 
satisfies,  12,at 

S12,at  =  -  <u  +  V)S12,aa-  (Ua  +  arv  )S12,a 


2  (S11  "  S22j]  uaa  ”  2(Sll,a 


-  S22,a)ua  +  1/2<S11  “  S22  +  2^)Va< 
+  1/2(t'll,a  -  S22,a)va  "  0r[vaS12,| 


(4 ,b. 14) 


+  vS12,aB  +  1/2(EH  '  S22  -  2^)UaB 

*  -  W*s] 
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while  for  B"0  (4.b.l4)  becomes 


S12,at  -  0 


(4.b.l4a) 


The  cross  derivative  of  the  deviatoric  stress  S,,  for  B>0, 
satisfies,  a, at 


aa 


S22,at  "  "  (u  +  arv)S22,aa  '  (ua  +  arva)S22,a 

-  +  arS12)uaa  -  arS12(aua  +  (S12  +  arfi)vQ 

(4.b.l5) 

+  S12,ava  "  Sr  [vaS22,B  +  vS22,aB  +  S12,auB 
+  S12uaB  “  3^  Vcb]  “  ^  “r 


while  for  B*0  (4.b.l5)  becomes 
S 


'22, at 


-  uS22,aa  ‘  uaS22,«  ‘  '  *r\,)  <4-b-15*> 


Now  we  state  the  results  for  differentiation  of  System  (4.b.2)- 
(4.b.8)  and  (4. b. 2a)- (4. b. 8a)  with  respect  to  B. 

The  continuity  equation  yields  for  the  cross  derivative  of  the 
density  pgfc,  for  B>0, 


Bt 


(vi  +  arv)paB  -  (us  +  arvB)po  -  P(J(ua  +  a^) 

-  p(uaB  +  °rVaB)  ‘  0r<2pBVB  +  PVBB  +  VPB8> 


-  (6rVpVB  +  VpB)  _  r  (pBV  +  pVB  “  "r  pv) 

-  (ar)6(vpa  +  Pva) 


(4.b.l6) 


while  on  the  axis  of  symmetry,  B=0,  (4.b.l6)  reduces  to 
p8t  =  -  2(6r>B  pV3 


(4.b.l6a) 


-37- 


Equation  (4.b.l6a)  is  not  used  since  (8r)g"0. 

The  cross  derivative  of  the  axial  velocity  ugt,  for  0>O,  satisfies 
u0t  “  -  (u  +  “rv)uo.0  "  (u0  +  arVua  +  ?  [Tll,a0 

Pg  “i 

+  “rS12,a0  -  p“  (Tll,a  +  arS12,o^J 
-  (v,  *  vu#s  -  i 2Ji  *  'Jagi)  ,4.b.r 


^  *  <S>,  («.  -  ^ 


(4.b.l7) 


for  0>O  while  for  0*0  (4.b.l7)  becomes 


u  _  2ffl  1  12'g 

0t  2  0  p 


(4.b.l7a) 


which  is  not  used  since  u,  as  well  as  p,  is  an  even  function. 

The  cross  derivative  of  the  radial  velocity  v^t,  for  0>O,  satisfies 

vflt  -  -  1“  +  <V',vt,B  ‘  |UB  +  Wv<. 

+  p  [[S12,<.8  +  “rT22,a6  P  IS12,a  *  “rT22,oi] 

-  Br  (-e-e  *  vvee  -  ^  )  ,4-b'191 


-  ‘Br>e(vvB  -  *  ? 


2S-0  „+S, 


(2S22+S11)  P0  (2S22+Sll)rf 


22  11' 
pr 
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while  for  @>0  (4.b.l8)  becomes 


*0t  "  uvo0  +  r(S12,oe  +  <V(iT22,a> 


"  0r(v6ve  ' 


*  2S22.BB  * 


(4.b.l8a) 


The  cross  derivative  of  the  internal  energy  eft.  ,  0>O,  satisfies 


the  tedious  relation 


B0t  "  '  <V  arv^ea0  '  <u0  +  arVea 


+  p  j(Tll,0  +  “rS12,0)ua  +  <T11  +  arS12,ua0 

+  (S12,0  +  arT22,0)va  +  (S12  +  arT22)vct0 

“  ^  [(T11  +  arS12)ua  +  <S12  +  arT22)v^J 

Pr  [v0e0  +  ve00  “  p(x22,0v0  +  x22v00  (4.b.l9) 

+  S12,0U0  +  S12U00  "  p^(x22v0  +  S12u0))] 

(6r)e  [vee  "  p(T22ve  +  si2ue)] 
pr  [VB  (p  +  S11  +  S22}  +  v(p0  +  Sllr0  +  S22,0) 

"  ^v<p  +  S11  +  S22>  "  v  <p  +  S11  +  S22)] 

"  (ar)6[Vea"  ^(S12Ua+  x22va>] 


which  simplifies  for  0=0 

x22Vfl 

egt  "  2'Br>e  0 


(4.b.l9a) 


Again  Equation  (4„b.l9a)  is  not  used  since  e  is  an  even  function 
around  0=0. 
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Sl2,Bt  =  -  (u  +  “rv)S12,a6  '  (UB  +  arVB)Sl2,a 

+  "7  (2p  “  S11  +  S22)ua0  "  “§  (S11,0  "  S22,0)ua 
+  1/2  (2y  +  Sl;L  -  S22)vag  +  1/2  (S11>g  -  s22,0)va 

"  er  [VSS12,g  +  vS12 , 03  ‘  1/2  (2y  “  S11  +  S22)u60 

*  V*  <S11,6  -  S22,6,Ub]  '  IBr’e  (VS12,B  <4'b-211 

-  1/2  (2„  -  SU  +  S22)ub)  -  «.r.6[vS12,a 

-  1/2  (2„  -  Su  +  S22)uo] 
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while  for  0*0  (4.b.21)  becomes 


S 


12,  et 


-  “s12,aB  +  1,1  +  -JVi>va8 

®ii  “®*>2 

+  (Vj1"  -  ^V^,u 


-8r[, 


V0S12,0 


a 

S, ,-S. 

T 


(4.b.21a) 


Finally,  the  cross  derivative  of  the  deviatoric  stress  S_2  g. 
for  0>O,  satisfies  ' 

S22,0t  *  "  (u  +  Qrv,S22,a0  "  (u0  +  arv0)S22,a 

,2U 


+  vS. 


_  (T“  +  arS12,ua0  "  arS12,Bua  +  (S12 

ar  rj  voB  +  S12,0vct  ■  Br  (V0S22,0 

4P  \ 

22, $6  +  S12,0U0  +  S12u00  "  T~  v00j 
-  (ve  -  rv>  -  <er>6  (vS22, 6 

+  S12u8  "  ve)  ~  (arJB  (vS22,a 
4y 


(4 .b. 22) 


+  S12ua  ‘  T  va) 


while  for  0=0  (4.b.22)  becomes 
S22,0t  =  (Br,0  T-  V0 


(4.b. 22a) 


Again  (4.b.22a)  is  not  used  since  S22  is  an  even  function  about  0=0. 

It  is  clear  then  that  Equation  (3.c»9)  may  now  be  generalized 
to  read 


wtt  "  P(3ct'  V  3aa  '  3a0'  300)w 


(4.b.23) 
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Since  w..  determined  from  Equation  (4.b.23)  is  multiplied  by 

2 

At  upon  substitution  into  Equation  <4.b.l)  it  is  only  necessary  to 
evaluate  the  finite  difference  approximation  of  the  spatial  differ¬ 
ential  operator  P  to  first  order  accuracy.  It  is  precisely  this 
fact  that  allows  us  to  use  uncentered  finite  difference  approxima¬ 
tions  to  spatial  derivatives  of  w  and  maintain  the  second  order 
accuracy  of  the  overall  scheme.  Conversely  it  is  precisely  the  lack 
of  a  fixed  regular  stencil  caused  by  the  boundary  of  the  domain 
cutting  the  Eulerian  mesh  which  leads  to  the  relaxation  of  requiring 
and  attaining  a  second  order  approximation  to  the  second  space  de¬ 
rivatives. 

Whenever  the.appropriate^neighbors  are  available  the  spatial 
derivatives  3a,  3^a,3g  and  3gg  are  approximated  by  centered  differ¬ 
ences  which  are  second  order  accurate.  For  example,  the  component 
u  is  approximated  by 


and  u  is  approximated  by 

(St>) 1 

and  the  error  made  through  this  approximation  is  of  the  order  of 
(Aa) 3 . 

Whenever  the  points  (i+1, j)  or  (i-l,j)  are  missing,  these 
formulas  are  replaced  by  noncentered  formulas  of  second  order  accur¬ 
acy.  For  example,  if  the  point  (i+l,j)  is  missing  and  the  distance 
between  the  point  (i,j)  and  the  boundary  (along  an  a  coordinate)  is 
a1>0  then  ua  is  approximated  by 


A(X  UVA 

ua  T  UB  +  oT[Ka-  ui,j  "  ui-l,j  (4>b'24) 


while  uaa  is  approximated  by 


u  ^2  < 
aa 


(a-^+Aa) 


_  Ui>  3 
a^Aa 


i~l» 


Aa  (a^+Aa) 


(4 .b. 25) 


In  these  formulas,  u0  is  the  value  of  u  at  the  boundary  point  B; 

it  is  defined  as  the  intersection  of  the  boundary  curve  with  the 
line  B=constant. 
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f 


?'■  '■ 


2 

Finally  £or  derivatives  cen  jring  is  not  possible;  how¬ 
ever  a  first  order  approximation  is  sufficient.  For  example 
uQ^  at  point  (a^,8j)  could  be  approximated  by 


i,j  :  ui-i>j  -  Ui,j-1  ;  ui-itj-i 


-  U. 


To 


~W 


This  formula  assumes  that  the  points  (i,j),  ( i— 1 ,  j )  and  (i-l,j-l) 
are  the  only  interior  points  available.  However,  it  is  not  necessary 
to  determine  all  the  interior  neighboring  points. 

In  order  to  simplify  the  logic  required  in  determing,  at  point 
(i,j),  which  interior  neighboring  points  are  available  for  use  in 
approximating  derivatives  by  difference  quotients,  it  is  convenient 
to  partition  the  search  procedure  into  two  phases:  a  nearest  neigh¬ 
bor  search  for  the  four  nearest  neighbors  followed  by  a  search  for 
the  four  remaining  neighboring  points.  Once  one  of  the  latter  points 
is  found  it  is  used  to  find  the  approximation  to  u  „  in  the  double 
Taylor  expansion  for  u  given  by 


u (a  +  Aa,  0  +  AB)  =  u(a,B)  +  Aau  +  ABu„ 

a  p 


+  I(A“2  uaa  +  2&aA(3  uaB  +  ^Sb* 


(4 ,b. 26) 


During  the  first  phase  of  the  boundary  search  about  the  point  where 
u(a,B)  is  interior  to  the  domain  of  integration  ua,  Uo,  u  and  Uog 
are  computed.  Once  it  is  determined  that  u(a+Aa,  +AB)  is  an  interior 
point  (here  Aa, AB  can  be  either  positive  or  negative) 
puted  from  Equation  (4.b.26). 

c)  Too  Near  Points 


u  „  is  com- 
aB 


There  are  a  set  of  interior  points  for  which  neither  the  two 
step  nor  the  one  step  difference  operators  can  stably  produce  an 
updated  solution  for  the  vector  w.  All  such  points  lie  within  a 
thin  annular  region  whose  boundaries  are  the  boundary  of  the  domain 
and  a  boundary  essentially  parallel  to  the  boundary  of  the  domain 
a  distance  on  the  order  of  one  third  Aa  or  AB.  All  points  lying 
within  this  band  take  values  interpolated  between  interior  and 
boundary  data. 


I 


J 


V 


Treatment  of  the  Boundary 


a)  Lagrangian  Representation  of  the  Boundary 

The  boundary  of  each  domain  is  considered  to  be  represented 
by  a  polygon  whose  vertices  are  Lagrangian  points  moving  with 
the  local  material  velocity.  At  each  point  i  the  differential 
equations 


dz.  dr. 

3tT  “  ui  >  3tT  =  vi  1  3  1'2"**»n  (5.a.l) 


are  solved  for  the  coordinate  pair  (z^,r^)  given  the  velocity 

vector  (^i)  at  the  point  i„  In  computational  space,  the  a-3  plane, 
i 

Equation  (5.a.l)  becomes 


da. 

cTt- 


u .  a 
l  z 


+  v.a 
i  r 


l  1,2, ... ,n 


(5.a.l*) 


df3i 

cTT 


ui6z  +  viBr 


The  values  of  velocity  at  each  point  on  the  boundary  is  known 
at  the  initial  time  so  that  the  new  boundary  position  is  computed 
to  be  the  polygon  with  vertices 


ai(t+At)  =  ai(t)  +  At(ui(t)az  +  vi(t)ar) 


(5.a.2) 


ei(t+&t)  =  ei(t)  +  At (Ui ( t)  sz  +  vi(t)3r) 

The  updated  values  of  the  velocity  components  at  the  new 
boundary  position  is  obtained  by  a  space-time  extrapolation  from 
the  interior.  The  data  chosen  for  extrapolation  to  the  kth  tracer 
particle  is  the  nearest  interior  neighbor  of  the  kt*1  point  with 
coordinates  (ak(t+At),  8k(t+At)),  i.eu 

uk(ak(t+At),  6k(t+At))  =  ufa^S-j^) 
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»*  if  i  tan  fri  i  fcj  i 


I"' 


vV» 


where  and  Bj  are  chosen  so  that 

(ak(t+At)  -  ai(t))2  +  (  6k(t+At)  -  B^t))2 

is  a  minimum  over  all  i  and  j ;  the  values  of  i  and  j  are  mesh 
crossings  interior  to  the  domain  of  interest. 

It  is  possible  to  correct  the  boundary  position  using  the 
latest  values  of  the  velocity  components  by  using  the  corrector 
formula 


ai(t+At)  =  ai(t) 


At 

T 


'X/ 

(u  (t) 


+  u (t+At) ) 


Bi (t+At)  =  0±  (t) 


+  ^  ( v  (t)  +  v  (t+At) ) 


(5. a. 3) 


'Xj  'b 

where  u  and  v  are  the  transformed  velocity  components  in  the  a-3 
plane;  they  are  given  by  the  right  hand  side  of  Equation  (5.a.l')> 

It  has  been  found  that  in  the  numerical  experiments  carried  out  thus 
far,  very  little  difference  appears  in  the  solution  of  the  position 
of  the  boundary.  This  is  probably  due  to  the  fact  that  the  time 
step  is  very  small  being  based  upon  the  relatively  high  sound  speed 
found  in  elastic  materials.  On  the  other  hand,  boundary  velocities 
are  usually  much  smaller  than  the  characteristic  sound  speeds  so 
that,  since  the  nearest  neighbor  is  usually  unchanged  for  the  kth 
boundary  point,  the  formula  (5. a. 2)  is  sufficient. 

The  boundary  is  moved  by  the  integration  of  the  differential 
equations  (5.a.l‘),  the  integrands  being  obtained  by  extrapolation 
from  interior  data.  As  this  boundary  sweeps  through  the  Eulerian 
mesh  data  must  be  defined  at  points  on  the  boundary  which  coincide 
with  the  Eulerian  mesh  lines.  This  set  of  points  is  used  to 
augment  the  set  of  interior  points  when  difference  approximations 
to  partial  derivatives  in  the  alpha  and  beta  direction  are  computed 
(see  Equations  (4.b.24)  and  (4.b.25)). 

b)  Free  Surface  Boundary  Conditions 

The  boundary  of  any  domain  is  composed  of  segments  which  con¬ 
stitute  one  side  of  a  slip  line,  i.e.,  an  interface,  or  a  free 
surface.  In  the  later  case,  a  purely  hydrodynamic  code  only 
requires  the  vanishing  of  the  pressure, 

p  =  0 

as  a  boundary  condition.  However,  for  elastic  domains,  the  pressure 
is  just  one  component  of  the  stress.  The  proper  condition  is  that 
the  normal  stress  must  vanish; 

Tn  =  0  (S.b.l) 


i 

I 


j 


! 

! 

i 
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A 


In  order  to  compute  Tn  we  let  tanijj  measure  the  local  slope  of 

the  boundary.  Let  fi  represent  the  local  unit  normal  vector  to  the 
boundary 


(siity 
COBlp 

and  €  represent  the  local  unit  tangent  to  the  boundary 


(5.b.2) 


*  /  -COSlJj  \ 

t  =  ( 5 ,b. 3 ) 

\  sintj/  j 


Let  the  stress  matrix  T  be  given  by 


/  T11  t12  ' 
\t12  x22 


Then  the  normal  stress  in  the  normal  dir'  ction  is 

J 

Tn*n  =  t11  sin  ip  +  2t12  sin^  cosij> 

2 

+  T22  cos  <l>  -  0 


( 5 .  b.  4 ) 


(5.b. 5) 


nere  Tij  =  -  p6iy 


Equation  (5.b.5)»  for  a  hydrodynamic  material  reduces  to  p=0. 
The  tangential  stress  in  the  tangential  direction 

Tt •  t  =  t: sin2i|j  -  2t12  sing»  cosip  +  t22  cos2ip  (5.b.6) 


may  be  arbitrary. 

In  addition  the  tangential  component  of  the  normal  stress 
vanishes 

/s  2  2 

Tn-t  =  (t22  “  x^^Jsinij;  cosg/  -  -r^2  (cos'ty  -  sin^) 

=  0  (5.b.7) 


-46- 


Under  conditions  (5.b.5),  (5.b.6)  ana  (5.b.7)  we  find  that 

A  A  A 

■  Tt*t  COS  Ip 

t12  ■  "  Tt*  t  sinip  cosip  (5.b.  8) 

A  A  m 

t22  ■  Tt* t  sin^tp 


are  the  desired  values  of  the  components  of  the  stress  tensor  in 
terms  of  the  extrapolated  values  of  the  tensor  T  obtained  from 
the  interior  of  the  domain. 

c)  Interface  Boundary  Conditions 

The  conditions  to  be  applied  to  the  boundary  at  the  interface 
are  simple  generalizations  of  the  three  conditions  described  in 
the  above  section. 

First  there  is  a  kinematic  condition  which  states  that  the 
jump  in  the  normal  velocity  at  the  interface  vanishes,  i.e. 

(a*1*  -  u*2hsinip  +  (v  ^  -  v^)cosip  =  0  (5.C.1) 


Here  we  have  used  superscripts  to  denote  the  material  number  on  each 
side  of  the  interface.  Condition  (5.^.1)  is  just  the  algebraic 
counterpart  of  the  jump  condition  [u*n]  =0  prescribed  at  the  inter¬ 
face. 


A  condition  on  the  stress  at  the  interface,  is  that  [Tn*n]  =  0 


i.e. , 


e(D+  cd) 

,22  bll 


-  P 


(1) 


e  (  2 )  ,  „(2) 

b22  fall 


+  ( 


S(1)-  S 
22  bll 


(1) 


S(2)-  S 
b22  bll 


(2) 


+P(2)) 


2  2 
-)  (cos  ip  -  sin  ip) 


(5.C.2) 


+ 


2(S{1)  - 


simp  cosip  =  0 


A  A 

In  addition  on  each  side  of  the  interface  Tn*t  =  0,  i.e.. 


(5.C.3) 


(Sj^”  sii^)  siniji  C0BI^  “  ^  (cos24<  -  sin2ii») 


^22^  ~  Sll^  sin^  COS^  ”  si2^cos2^  ”  sin2^) 


=  0  (5.C.4) 


The  quantity  Tn*n  is  computed  from  each  side  of  the  interface. 
In  order  to  assure  that  the  normal  velocity  remains  continuous  we 
introduce  a  new  value  of  Tif»n  in  terms  of  these  computed  value  via 


Tfi-rf 


(1) 


Tn*n 


(2) 


“nr 


(2)  a(1) 

o  Tn-n 

T2V  - 


+  P 


(5.C.5) 


A  value  of  Tn*t  is  formed  via  the  welded  boundary  condition 
assumption 

Tn*  t (1)  =  Tn.€<2)  =  1/2 (Tn-  £(1)  +  Tn-t(2))  (S.c.6) 

The  value  of  Tt*£  is  obtained  from  the  interior  of  each  side 
of  the  interface,  i.e.  both  Tt-t'1',  T£*t'2'  are  extrapolated 
along  almost  characteristic  directions  from  the  corresponding 
interior  of  each  material. 

We  now  have  the  fact  that  given  Tn»n  from  Equation  (5.C.5) 
and  Tn«t  from  the  welded  ^boundary  condition  (5.C.6)  as  well  as  the 
extrapolated  values  of  Tt»t  from  each  side  one  can  compute  the  stress 
components  from  the  set  of  linear  equations 

A  A  2  2 

Tn»n  +  p  =  2siniJ>  cosi}/  S12  +  sin  ^  +  cos  ^  S22 

Tn«tf  =  -  (cos  f  -  sin  i)j)3^  -  sin^  costy 

-  sin^  cosip  S22 

Tt*t  +  p  =  -  2sintj;  cosij)  S^2  +  cos2^  +  sin2ij,  S22 

The  determinent  of  the  above  system  is  unity. 
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By  applying  Kramer’s  rule  we  find 


S 

S 

S 


11 

12 

22 


A  A  2  A  A  2 

=  Tn*n  sin  ip  +  Tt*t  cos  <p 

-  2Tn*t  sintp 

o 

o 

00 

+  P 

(5.C.7) 

*  (sin^ip  -  cos^ip)  Tn?t  + 

(Tn?n  -  Tt«t) 

sintp 

CO  Sip 

(5.C.8) 

A  A  9  A  A 

A  A 

2 

=  Tn*n  cos  tp  +  2Tn?t  sinip 

COSIp  +  Tt‘t 

sin  tp 

+  p 

(5.C.9) 

The 
taken  to 


condition  on  the  normal  velocity  of  the  interface  u.  is 
be 


Uj_  =  <p(1)  a*1* 


P(2) 


P(2)) 


(5.c. lOi 


while  the  tangential  component  uH  ^  ,  u,/  2  are  extrapolated  from 
the  corresponding  interior  positions  along  almost  characteristic 
directions. 

d)  Characteristic  Equations 

We  have  carried  out  a  formulation  of  the  equations  of  motion 
at  interface  boundaries  in  characteristic  form.  It  is  desirable 
to  formulate  the  difference  problem  on  a  coordinate  system  in  which 
the  two  coordinates  are  locally  orthogonal  and  parallel  to  the 
boundary.  Since  the  boundary  will  exhibit  a  spatial  variation  in 
slope  which  will  change  in  time,  it  is  convenient  to  recast  the  basic 
differential  equations  in  a  form  where  differentiation  is  automatic¬ 
ally  carried  out  in  a  direction  more  or  less  normal  and  parallel  to 
the  boundary.  The  corresponding  difference  scheme  thus  generated 
can  be  aligned  with  the  boundary  -  differences  being  taken  perpendi¬ 
cular  and  parallel  to  the  boundary. 


One  first  takes  an  appropriate  linear  combination  of  the 
momentum  equations  and  constitutive  relations,  i.e.. 


P5?  +  P2  -  <S 


S12 

11, z  +  S12 ,  r  +  ~T~) 


T  Dv  .  „  ,  2S22+SllH 

LPDt  +  pr  (S12,z  S22,r  r 


-  tan9 

sin9 

/p/p 


/p/p  cosg 


(5.d. 1) 


sin9 

/p/p 


\Pl  .  211  (2Vr  .  Uz  .  S,]  . 
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Then  by  introducing  the  definition  of  characteristic  differentiation 
in  terms  of  the  particle  derivative  via 


1/2 


3 


/2 

+  COS0 


3 


3t'5t  +  (^8inu  j*  T  ^JCOBO  57 

,V2 


-  k  +  (u  +  (p)sine  >  k  +  (v  +  (^)c/°2se)  k 


(5.d  . 2) 


Equation  (5.d.l)  can  be  written  as 

-dS 


tdu  -  dv  . 
jjr  cosfl  -  tz  sinO 


3t 


sin20 


rdSll  dS22~*[  cos20  dS12 

Lat  dt  ,—TZ  dt 


^  jT7o  Ldt  dt  J  J~V7t 


slnf)  “  v  cos9)  +  ^  (sn  cos^B  +  Sj2  sin*1 


-  p  -  2  S12  sin0  cos0)  +  —  (S12  cos0  -  ( 2S22  +  S^)  sin0) 


(5.d. 3) 

This  is  the  characteristic  equaticr^along  the  conoid  with  the  local 
disturbances  propagate  with  speed  =  p/p  ,  the  shear  speed,  rela¬ 
tive  to  the  particle  speed.  The  angle  0  is  measured  from  the  image 
of  the  characteristic  on  the  plane  t=0  from  the  r  axis  and  the 
directional  derivative  d/d£,  defined  in  the  plane  t=0  and  directed 
tangent  to  the  base  of  the  conoid,  is  defined  by 


k  =  cose  k  ■  sin9  h 


(5.d. 4) 


In  addition  to  the  shear  conoid,  there  appears  another 
characteristic  cone  defined  from  the  characteristic  speed 
2  2 

C,  =  C  +  4p/3p,  C  being  the  sc  5  speed.  Obviously,  disturbances 
which  propagate  at  speed  Cd  relative  io  the  material  particle  speed 
will  travel  further  in  time  At  than  an  disturbances  which  travel 
at  the  shear  speed  Cg.  Hence,  for  wave  motion  in  an  elastic  medium, 
the  shear  conoid  lies  inside  the  sonic  conoid. 

In  order  to  find  the  characteristic  equations  along  the  sonic 
cone  we  again  take  linear  combinations  of  the  continuity  and  momentum 
equations,  compatability  equations  and  the  particle  path  equation 
which  relates  changes  in  pressure  to  density  on  a  particle.  The  re¬ 
sult  is 
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-  <’u..  *  s12,r  *  ^3  *  [-1  *  pr 

'  (s12,z  *  S22,r  ~  r  U>]  *  JT  [St  "  c2  if] 

-  [b^  -  3“‘**z  -  vr  -  ?>]  '5-0 

-  W-ivaa-  -  n<ur  ♦  .,,]  -  ss^i  [£a 


(5.d. 5) 


-  ¥■  (2vr  -  “z  '  ?>] 


Equation  (5.d.5)  can,  after  a  fair  amount  of  algebraic 
manipulation  is  performed,  be  put  into  the  characteristic  form 


cVi+s?  o 


£  sine  -  U  cose]  *  |f 


2  dSll  ^S12  2  dS22 

-  sin  9  g^ - 2  sine  cose  g^ - cos  e  gfc 

=  (j^  -  pC2)  g^-  [u  cose  -  v  sinej  (5.d.6) 

+  C  / 1  +  ^ j  (sine  cose  (St,-  S,0) 


'11  22 ' 


+  (cos2e  -  sin2e>  +  |(Cy  1  +  ~^~2  Jsj_2  sine 


+  (2  S22  +  SX1 


Jcosej  +  -  pC2)v>/r 


We  have  used  the  following  definitions  in  Equation  (5.d.6) 
The  derivative  in  the  bicharacteristic  direction  is 

at =  St  +  C71  +  ^2  lh  sine  +  h  cos9)  (5-d-7) 
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while  the  material  particle  derivative  is  given  by 


Dt“TE+ull+vlr  (5.d.C) 


As  before  Equation  (5.d.4)  defj  •  the  directional  derivative 
d/d(  which  is  used  in  Equation  (5.d.  l 

The  method  of  computation  is  straight  forward;  the  derivatives 
appearing  in  the  above  equations  are  replaced  by  differences  in 
the  direction  in  which  the  derivatives  were  defined.  For  example. 
Equation  (5.d.3)  can  be  written  in  finite  difference  form  as 


P (cos61  (1)  u*1)  -  Sine^1*  v*1*) 

-  2  sine^11  cos01  (1>  {S^1-  > 


(cos  0 


2«  (1)  _  ,4-2n  UK  eU) 


-  sin  0^  )  S|2 


(5.d.9) 


■J  P^V11  (coso/1)  „<»>  -  .me/1’  v(0)l 


-  2  sin6 


(1)  cosO^11  (S^*  -  S ^ 


-  (cos201(1)  -  sin201{1) 


+  At.  (difference  approximations  to  the  RHS  of 
Equation  (5.d.3)) 


The  differencing  is  along  a  generatrix  of  the  shear  cone  and 
along  a  line  perpendicular  to  the  generatrix  in  the  z-r  plane. 

This  point  of  intersection  is  denoted  by  the  superscript  zero  while 
the  vertex  of  the  shear  cone  is  denoted  ry  superscript  one.  The 
subscript  one  denotes  the  particular  angle  theta  chosen  to  define 
the  ray  for  the  integration;  the  superscript  on  the  angle  0  denotes 
the  material  number.  Equation  (5.d.9)  as  it  now  stands  has  been 
written  on  one  side  of  the  interface  and  therefore  a  similar  relation 
must  be  written  for  the  other  side.  Thus,  whenever  we  have  a  term 


(1) 


it  is  replaced  by  e, 

4  At  /It  * 


replaced  by  p  ^  and  u^ 


(2)  (1) 


A 


is  replaced  by 


(2) 

p0 


(1) 


IS 


and  v  are  understood  to  be  the  velocity 
components  at  the  vertex  of  the  conoid  on  the  other  side  of  the 
interface,  i.e.  u^2)  and  v'2).  This  gives  us  a  second  algebraic 
equation  representing  integration  along  the  shear  cone  to  the  inter¬ 
face  from  the  second  material. 
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The  lame  procedure  ii  now  applied  to  Equation  (5.d.6)  with  the 
exception  that  two  angles  along  the  dialatational  conoid  are  chosen. 
In  general  these  angles  are  distinct  from  those  chosen  for  the  shear 
conoid.  This  gives  us  four  more  algebraic  relations  for  the  twelve 


unknowns  u 

S<2) 


(1) 


11 


c  (2) 
S12  ' 


v(1>, 


j  (2) 

>22 


P(1), 


e  (1) 
S11  ' 


e(l) 
S12  ' 


5d) 

>22 


and  u<2>, 


v<2>, 


P(2), 


One  now  looks  at  the  basic  stress  strain  relations  which  include 
the  terms  accounting  for  rotation: 


S 


11,  t 


S 


12, t 


S22,t 


(u  S 


11,  z  +  S12  vz  ’  7  uuz> 


lv  sll,r  "  S12  ur  +  T1  vr>  - 
(u  S10  „  -  (U  +  vz 


'12, z 
512,r 


®ll“S22 

(v  S10  ^  -  <M  -  ■  ~~~2  >  V 


(u  S22,z  "  S12  vz  +  ¥  uz> 


(v  S22,r  +  S12  Ur  ’  3^  V  " 


2m  v 
3“  ? 


2m  v 
3  r 


(S.d.lOa) 


(S.d.lOb) 


(5.d. 10c) 


If  we  neglect  these  rotation  terms  and  multiply  Equation 
(5.d.l0a)  and  Equation  (5.d.l0c)  by  3/2  they  may  be  added  to 
obtain 


7  at  (S11  +  S22}  -  n<uz  +  vr  -  T-)  (S.d.ll 

Equation  (S.d.lOb)  may  now  be  added  and  subtracted  to  Equation 
(5.d.ll)  to  yield 


d 

3t 


<SU  *  S22>  ±  ht]  ■  “  [<fi  i  fr>  <“  i  »>*  t]  <5-d-12' 


Equation  (5.d.l2)  has  a  directional  derivative,  gp  defined 

along  the  particle  path  while  the  right  hand  side  is  expressed  in 
terms  of  a  fixed  directional  derivative,  i.e.  fixed  along  the 
directions  +  45  degrees.  Integration  of  Equation  (5.d.l2)  along 
the  particle  path  on  each  side  of  the  interface  gives  us  two  more 
independent  algebraic  relations  for  the  twelve  unknowns.  The 
remaining  four  relations  are  obtained  from  the  boundary  conditions. 
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If  the  slope  of  the  interface  is  tan*  ■  dz/dr  then  the  con¬ 
tinuity  of  normal  velocity  yields 


(u(1>  -  u<2>)  sin*  +  (v<1>  -  v<2>)  cos*  ■  0 


(5.d.l3) 


The  remaining  conditions  are  obtained  from  constraints  on 
the  stress  tensor  T: 


'll 


'12 


12 


22 


(5.d.l4) 


Let  n  = 


sin* 
cost}) . 


be  the  unit  normal  to  the  interface  with  slope 

1  (  A  (  — cos 

tan*;  therefore  the  unit  tangent  is  t  =  \  ,, 

a  ^  \  sir^ 

the  normal  stress,  [Tn*n]  =  0  yields  our  second  ' 

Equation  (5.C.2). 


Continuity  of 
boundary  condition. 


The  remaining  two  conditions  can  be  obtained  by  placing  a  con¬ 
dition  on  tjie  components  of  the  normal  stress,  ^n,A  in  the  tangential 
direction,  t.  We  again  invoked  the  condition  Tn*t  *  0  on  side  1 
and  side  2. 


Hence 

sin*  cos*  -  (cos2#  -  sin2*)  =  0  (5.d.l5) 

-(s|2)-  sj2))  sin*  cos*  -  sj2) (cos2*  -  sin2*)  =  0  (5.d.l6) 
are  the  required  conditions. 

The  entire  system  of  equations  can  be  written  in  the  following 

form 

A  w  =  b  (5.d.l7) 
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where  the  unknown  vector  w  has  entries  given  by 


while  the  coefficient  matrix  A  of  order  12  has  row  entries  defined 
by 

1)  Continuity  of  normal  velocity.  Equation  (5.d.l3) 

2)  Continuity  of  normal  stress.  Equation  (5.C.2) 

3)  Boundary  condition.  Equation  (5.d.l5) 

4)  Boundary  condition.  Equation  (5.d.l6) 

5)  Difference  equations,  Shear  conoid  material  1, 

Equation  (5.d.9) 

6)  Difference  equations,  Shear  conoid  material  2 

7)  Difference  equation,  Dialatational  conoid-material  1 
with  angle  G^d),  (difference  form  of  Equation  (5.d.6) 

8)  Difference  equation,  sonic  conoid-material  1 
with  angle  82'^ 

9)  Difference  equation,  sonic  conoid-material  2, 
with  angle 


10) 


Difference  equation,  sonic  conoid-material  2, 
with  angle  Gj'2' 
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11)  Difference  equation,  Particle  path  relationship  for 
material  1,  Equation  (5.d.l2)  in  difference  form 

12)  Difference  equation,  Particle  path  relationship  for 
material  2,  Equation  (5.d.l2)  in  difference  form 


In  the  above  listing  the  entries  actually  used  in  row  8  and 
row  10  are  the  differences  of  the  two  difference  equations  defined 

by  integration  along  the  angles  0^,  for  each  of  the 
materials,  i>l,2. 

The  above  formulation  has  been  programmed  and  is  presently 
undergoing  testing.  It  is  considered  to  be  a  prototype  of  methods 
which  may  be  used  to  predict  interactions  of  two  dissimilar  mater¬ 
ials  undergoing  elastic  impact. 
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VI. 


Results 


In  order  to  demonstrate  the  versatility  of  the  algorithm, 
two  problems  were  considered.  The  first  problem  consists  of  a 
penetrator  composed  of  a  90/25  tungsten  alloy  impacting  at 
0.142  cm/ y sec  upon  a  one  inch  thick  RHA  plate.  The  configuration 
at  moment  of  impact  is  shown  in  Figure  (1).  The  second  problem, 
shown  in  Figure  (2),  is  a  similar  tungsten  projectile  but  now 
enclosed  in  a  Maraging-300  steel  sheath.  The  impact  velocity  is 
the  same  as  in  the  first  problem.  Each  material  in  both  problems 
is  assumed  to  have  an  equation  of  state  which  is  qiven  by 
Tillotson: 

Compressed  states:  p*p 

c 

P_  =  (a  +  -=£ - )EP  +  AC  +  BC2 

C  —.+1 


where  C  =  n  -  1  ,  n  =  p/pq 


P  >  PQ 
0  <  E  <  E 


Expanded  States:  p=p 


aEP  ♦  /-§5£ 


+1 


+  ACe 


V 


1,  -«c± 


i)' 


o  <  p 


0 

E  >  E 


s 


Intermediate  States: 


P  = 


(E  -  Eg)pe  +  (e;  -  E)p„ 


E'  -  E 
s  s 


P  <  P, 


E  <  E  <  E' 
s  s 


The  constants  used  in  the  above  equation  of  state  for  the 
present  calculation  are  given  in  the  following  table. 


The  coefficient  of  viscosity  K  in  Equation  (4. a. 7)  was 
taken  to  be  1.8.  The  two  transformations  (2.1)  and  (2.2),  which 
are  presently  coded  inline  rather  than  in  a  function  subprogram! 
form,  were  used  in  each  material.  The  constants  for  the  trans¬ 
formations  applied  to  each  material  domain  are  given  in  the 
accompanying  table. 


L 

Transformation 

__rmax 

l^T72 

(2.1) 

d 

Transformation 

D 

(2.2) 

q 

90/25 

Tungsten 

17 

1.75/16.5 

1.75 

Az'1 

0 

RHA 

24 

9/23.5 

0.5 

Az"1 

0 

Maraging-300 

17 

1.75/16.5 

1.75 

Az"1 

0 

The  constants  used  in  the  transformations  for  the  projectile 

were  chosen  with  d>r  so  that  uniform  spacing  in  r  is  achieved. 

—  max 

For  the  target,  fine  soacing  in  the  region  of  impact  is  desired  so 
d<  r. 


max 


The  mesh  (iAa,  j  A3) ,  which  was  chosen  such  that 
1  <_j  <_J ,  is  given  in  the  following  table  for  each  of 


1  <i  <_I, 
the  domains. 


I 

J 

ig-right 

j0 

Vleft 

90/25 

Tungsten 

121 

20 

120 

17 

21 

RHA 

31 

25 

25 

24 

16 

Maraging-300 

61 

20 

60 

17 

11 

Here  iQ  and  jQ  represent  respectively  the  initial  a  position 

of  the  material  right  justified  on  the  mesh  and  the  initial  maximum 
height  of  the  material  in  the  6  direction.  The  initial  starting 
value  of  ig  left  justified  is  also  shown.  The  minimum  value  of  jg 

is  1  for  all  domains. 
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The  final  configuration  after  30  microseconds  for  the  un¬ 
sheathed  penetrator  is  shown  in  Figure  3.  The  letter  P  appears 
at  each  mesh  point  when  the  material  is  at  a  stress  level  which 
satisfies  the  yield  condition.  At  such  points  the  material  be¬ 
haves  plastically.  Figure  4  corresponds  to  the  final  state  of 
the  sheathed  penetrator  at  30  microseconds.  At  this  time  there 
is  slightly  less  deformation  of  the  nose  of  this  penetrator  com¬ 
pared  to  the  unsheathed  penetrator.  Because  of  the  velocity  and 
density  ratio  between  the  penetrator  and  target, the  rod  is  not 
significantly  deformed  during  the  penetration  process  except  at 
the  nose  where  some  blunting  and  mushrooming  of  the  initial  spheri¬ 
cal  shape  takes  place. 

In  both  cases/  Figures (3)  and  (4)  we  see  that  the  rod  has 
completely  penetrated  the  target  and  is  emerging  from  the  far  face. 
Since  spall  and  other  fracture  mechanisms  have  not  been  incorpor¬ 
ated  into  the  present  model,  plugging  is  not  accounted  for;  ob¬ 
viously  the  back  face  of  the  target  has  stretched  beyond  that  which 
would  occur  for  RHA. 

Figure  (5)  compares  the  pressure  history  at  the  interface  on 
the  axis  of  symmetry  of  the  penetrator  for  both  of  the  above  prob¬ 
lems,  while  Figure  (6)  compares  the  pressure  history  in  the  target 
for  both  problems.  The  □  symbol  is  the  two  material  curve  while  the 
A  is  the  three  material  curve.  For  this  high  velocity,  the 
pressure  is  continuous  across  the  interface  since  the  stress 
revel  is  on  the  order  of  1/10  the  pressure  levels.  Thus,  on  the 
scale  being  plotted,  both  sets  of  curves,  penetrator  and  target, 
are  nearly  coincident.  Figure  (7)  is  a  plot  of  the  penetration 
depth  D  vs  time.  Over  the  first  20y  seconds  the  penetrator  leading 
edge  is  moving  at  an  average  speed  of  -.06  cm/ysec.  It  is  clear, 
from  this  graph,  that  the  speed  is  increasing  since  the  target  is 
failing  beyond  this  time.  Similar  behavior  is  exhibited  for  both 
problems.  At  30y  sec.,  the  sheath  is  still  moving  at  -0.142  cm/psec. 
as  is  the  penetrator  for  approximately  the  last  90%  of  its  length. 
Hence  the  residual  velocity,  although  not  computed  by  integration 
over  the  volume  region  of  space  defining  the  projectile  is  on  the 
order  of  >.9  initial  velocity. 

Computed  stress  levels  in  the  sheath  at  this  time  are  low. 

The  maximum  pressure,  at  the  interface  near  the  leading  edge,  is 
approximately  one  kilobar  but  on  the  average  the  pressure  lies  be¬ 
tween  one  kilobar  and  one  hundred  bars.  Maximum  and  average  stress 
levels  are  similar  The  maximum  pressure  transmitted  to  the  target 
is  somewhat  larger  for  the  sheathed  rod  impacting  although  the 
pressure  histories  are  similar.  In  both  cases,  the  projectile  re¬ 
mains  compressive  near  the  leading  edge  out  to  25  microseconds. 

For  the  unsheathed  penetrator  problem,  approximately  250  mesh 
points  were  used  in  the  targets  and  2000  points  in  the  penetrator. 
This  problem  ran  to  30  microseconds  problem  time  in  2271  seconds  on 
a  CDC  6600.  For  the  problem  with  sheath  the  target  and  the  pene¬ 
trator  have  the  same  number  of  points  as  the  first  problem,  while 
the  sheath  contained  1000  mesh  points.  The  computation  time  for 
this  problem  is  2169  seconds. 
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